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SUMMARY 

This  is  a  four-part  final  report  on  the  research  supported  by  the  Air  Force  Office  of 
Scientific  Research  Center  under  Grant  F49620-99- 1-0037,  entitled  Real  Time  Predic¬ 
tive  Flutter  Analysis  and  Continuous  Parameter  Identification  of  Accelerating 
Aircraft. 

1.  Motivations  and  research  objectives 

Flutter  clearance,  which  is  part  of  any  new  aircraft  or  fighter  weapon  system  de¬ 
velopment,  is  a  lengthy  and  tedious  process  from  both  computational  and  flight  testing 
viewpoints.  An  automated  approach  to  flutter  clerance  that  increases  flight  safety  and 
^educes  flight  hours  requires  as  a  stepping  stone  the  development  of  a  real  time  flutter  pre¬ 
diction  capability.  Such  a  fast  analysis  tool  can  be  designed  if  the  coupled  fluid/structure 
aeroelastic  system  is  represented  by  a  simplified  mathematical  model  that  can  be  quickly 
adapted  to  changes  in  flight  atmospheric  conditions,  aircraft  mass  distribution  (weapon 
systems),  fuel  loading,  and  Mach  number,  and  if  the  current  parallel  processing  technology 
is  exploited. 

Furthermore,  flight  testing  is  always  required  to  establish  the  flutter  envelope  of  an 
aircraft.  The  traditional  method  for  determining  such  an  envelop  uses  test  data  extracted 
from  the  vibration  response  of  the  aircraft  at  fixed  flight  conditions.  The  aircraft  is  first 
trimmed  to  a  specific  flight  condition  (Mach  number  and  dynamic  pressure),  then  its 
aeroelastic  response  is  deliberately  excited  by  applying  an  input  to  a  flight  control  surface. 
The  frequency  and  damping  of  the  excited  aeroelastic  response  are  typically  extracted  from 
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the  vibration  data.  By  repeating  this  test  at  many  flight  conditions,  the  flutter  envelope 
can  be  determined.  Such  a  traditional  approach  requires  that  the  aeroelastic  response 
be  measured  at  many  different  flight  conditions.  This  often  requires  a  large  number  of 
flight  test  hours,  a  process  which  not  only  costs  money  but  also  exposes  test  pilots  to 
proportionately  increased  risk.  However,  this  test  procedure  can  be  expedited  if  data 
collected  from  continuously  varying  flight  conditions  can  be  used  to  extract  the  needed 
flutter  Humping  and  frequency  values  from  an  accelerating  flight  profile.  In  that  case,  it 
may  be  possible  to  greatly  reduce  the  number  of  flight  hours  required  for  establishing  the 
flutter  envelope. 

The  Air  Force  Flight  Test  Center  at  the  Edwards  Air  Force  Base  (AFB)  has  expressed 
great  interest  in  the  above  two  problems,  and  therefore  we  have  proposed  to  conduct 
a  three-year  research  effort  in  real  time  flutter  analysis,  and  the  continuous  parameter 
identification  of  an  accelerating  aircraft.  More  specifically,  we  have  proposed  to  develop  a 
simplified  flutter  analysis  method  that  can  be  run  real  time  to  provide  predictive  frequency 
and  damping  values  for  maneuvers  as  flown.  The  enabling  technology  of  such  a  real  time 
flutter  analysis  capability  is  a  formulation  of  the  aeroelastic  problem  that  allows,  among 
other  things,  partial  pre-solutions  and  the  usage  of  parallel  processing.  We  have  reported 
on  this  research  activity  dining  the  first  year  of  support. 

We  have  also  proposed  to  develop  a  parameter  identification  technique  that  can  be  used 
to  extract  frequency  and  damping  values  of  an  aircraft  that  is  continuously  accelerating. 
This  technique  should  reduce  both  the  cost  and  risk  involved  in  determining  the  flutter 
envelopes  of  fighters.  Here,  we  report  on  this  effort  which  has  been  conducted  during  the 
last  two  years  in  collaboration  with  the  researchers  and  engineers  of  the  Air  Force  Flight 
Test  Center  at  the  Edwards  AFB. 

2.  Major  accomplishments  during  the  last  two  years  of  support 

We  have  determined  that  two  fundamental  issues  must  be  addressed  before  a  method 
for  the  continuous  parametric  identification  of  an  accelerating  aircraft  can  be  developed. 
The  first  issue  deals  with  how  the  aeroelastic  properties  of  an  aircraft  can  be  affected  by 
a  constant  acceleratic  j.  in  a  level  flight  or  during  maneuvering.  In  particular,  is  it  possible 
to  relate  in  a  simple  way  the  aeroelastic  parameters  measured  in  an  accelerated  flight  to 
those  measured  in  stabilized  flight  conditions?  The  second  issue  is  related  to  the  fact  that 
most  if  not  all  identification  methods  used  in  practice  implicitly  assume  that  the  given 
aeroelastic  system  is  linear  and  non-varying  in  time.  Whether  these  methods  can  still 
be  used  to  analyze  accelerated  flight  data,  or  whether  new  methods  are  required  for  this 
purpose  was  an  open  question. 

During  the  last  two  years  of  funding,  we  have  addressed  important  aspects  of  the 
above  two  issues  by  performing  appropriate  analytical  studies  and  CFD  based  numerical 
simulations.  More  specifically,  we  have  first  considered  a  typical  NACA  0012  wing  sec¬ 
tion,  and  investigated  the  effects  of  horizontal  and  vertical  accelerations  on  the  aeroelastic 
response  of  this  system.  We  have  shown  analytically  that  accelerating  the  wing  section 
introduces  additional  inertia  forces  and  modifies  the  torsional  stiffness  of  the  aeroelastic 
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system  by  a  negligible  quantity.  Next,  we  have  developed  a  procedure  for  extracting  the 
frequencies  and  damping  coefficients  of  an  aeroelastic  system  from  time  history  data  of  an 
accelerated  flight  profile.  This  procedure  is  based  on  a  modification  of  the  ERA  algorithm, 
and  its  interaction  with  a  windowing  concept  where  the  Mach  number  is  assumed  to  be 
constant  in  each  window.  In  order  to  validate  this  new  procedure,  we  have  upgraded  our 
computational  aeroelasticity  capability  to  handle  accelerated  flight,  which  was  by  itself 
an  interesting  and  rewarding  research.  We  have  simulated  both  cases  of  stabilized  flight 
conditions  and  accelerated  flight.  We  have  compared  the  generated  results  and  formulated 
conclusions  regarding  the  theoretical  and  practical  feasibilities  of  extracting  the  flutter 
envelope  of  an  aircraft  from  an  accelerated  flight  data.  These  conclusions  essentially  high¬ 
lighted  the  significant  potential  of  this  new  concept  of  flight  testing.  Motivated  by  our 
success  for  the  NACA0012  airfoil,  we  have  repeated  our  simulations  of  the  continuous  pa¬ 
rameter  identification  of  an  accelerating  aeroelastic  system  for  a  typical  F-16  wing  section. 
We  have  designed  this  wing  section  from  geometrical  and  structural  data  obtained  from 
the  Edwards  Air  Force  Base.  We  have  simulated  the  continuous  parameter  identification 
for  the  F-16  Block  40  typical  wing  section  in  accelerated  flights  with  up  to  0.05  Mach 
per  second  and  for  flight  regimes  extending  from  subsonic  to  supersonic.  We  have  shown 
that  the  aeroelastic  parameters  identified  in  accelerated  level  flight  at  a  given  altitude  are 
almost  identical  to  those  obtained  in  stabilized  flight  conditions,  which  justifies  the  pro¬ 
posed  new  concept  of  flight  testing.  However,  we  could  not  match  perfectly  the  results  of 
the  numerical  simulations  using  the  typical  wing  section  with  those  of  the  actual  test  of 
the  F-16  Block  40,  particularly  for  damping.  This  was  because  the  available  test  data  is 
for  a  loaded  wing,  whereas  our  typical  wing  model  was  derived  from  a  clean  wing  model, 
and  because  the  typical  wing  section  approach  is  valid  only  for  uniform,  straight  and  high 
aspect  ratio  wings.  All  these  developments  and  finding  are  described  in  the  attached  paper 
“Expanding  a  Flutter  Envelope  Using  Accelerated  Flight  Data:  Application  to  an  F-16 
Fighter  Configuration”. 

Next,  we  have  designed  a  three-dimensional  F-16  Block  40  aeroelastic  model  using 
two  incomplete  finite  element  structural  models  acquired  from  Lockheed-Martin:  (1)  a 
static  model  which  does  not  contain  the  mass  distribution,  and  (2)  a  linear  “fish-stick” 
dynamic  model  which  contains  the  needed  mass  information  but  which  is  not  adequate  for 
stress  analysis.  Using  this  model,  we  have  obtained  numerical  results  for  the  F-16  Block 
40  that  agreed  amazingly  well  the  flight  test  data.  This  latest  development  is  described 
in  details-  in  the  attached  paper  “‘Application  of  a  Three-Field  Nonlinear  Fluid-Structure 
Formulation  to  the  Prediction  of  the  Aeroelastic  Parameters  of  an  F-16  Fighter”. 

3.  Impact  on  the  state-of-the-art  of  flight  testing 

Motivated  by  the  results  we  have  generated  under  this  Grant,  the  Air  Force  Test 
Pilot  School  performed  on  May  15-16,  2000,  two  suites  of  accelerated  flight  tests  designed 
by  the  Principal  Investigator  and  his  collaborators  at  Edwards.  The  first  series  of  tests 
was  for  a  clean  (no  stores)  F-16  configuration,  and  the  second  for  an  F-16  configuration 
known  to  cause  Limit  Cycle  Oscillations  (LCO).  In  both  cases,  the  flight  test  team  was 
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able  to  accomplish  the  set  objectives  using  2.5  times  less  sorties  than  when  using  the 

classical  stabilized  flight  testing  approach,  which  demonstrates  the  potential  impact  of  our 

accomplishments  on  the  flight  testing  technology. 
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Expanding  a  Flutter  Envelope  Using 
Accelerated  Flight  Data:  Application  to  an 
F-16  Fighter  Configuration 


f  C.  Farhat  *  f  C.  Harris,  and  \  D.  J.  Rixen  ■ 
f  Department  of  Aerospace  Engineering  Sciences 
and  Center  for  Aerospace  Structures 
University  of  Colorado,  Boulder,  CO  80309-0429 
t  LTAS-Dynamics  of  Structures 
University  of  Liege,  Belgium 

Using  both  simple  analytical  investigations  and  complex  numerical  simulations,  we 
show  that  typical  accelerations  of  an  aircraft  do  not  affect  significantly  its  aeroelastic 
parameters.  This  suggests  that  flutter  testing  could  be  performed  in  accelerated  flight, 
thereby  reducing  the  cost  and  risk  involved  in  determining  the  flutter  envelopes  of  fighters. 
We  illustrate  this  concept  with  numerical  simulations  for  an  F-16  fighter  configuration 
and  compare  the  obtained  results  with  flight  test  data. 


1  Introduction 

Flight  testing  is  always  required  for  establishing 
the  flutter  boundary  of  an  aircraft.  The  traditional 
method  for  determining  a  flutter  envelope  uses  data 
extracted  from  the  vibration  response  of  the  aircraft 
at  stabilized  flight  conditions.  Typically,  the  aircraft  is 
first  trimmed  to  a  specific  Mach  number  and  dynamic 
pressure,  then  its  aeroelastic  response  is  excited  by 
applying  an  input  to  one  of  its  control  surfaces.  The 
frequency  and  damping  of  this  response  are  then 
extracted  from  the  vibration  data.  By  repeating  this 
test  at  several  different  stabilized  flight  conditions, 
the  flutter  envelope  is  determined. 

The  traditional  testing  approach  summarized  above 
implies  a  relatively  large  number  of  flight  test  hours,  a 
process  which  is  not  only  expensive,  but  also  exposes 
test  pilots  to  increased  risks.  One  way  to  expedite  this 
test  procedure  is  to  develop  a  method  for  expanding 
the  flutter  envelope  of  an  aircraft  which  can  use  data 
collected  from  continuously  varying  flight  conditions. 
By  extracting  the  critical  damping  and  frequency 
values  from  an  accelerating  flight  profile,  it  may 
be  possible  to  substantially  reduce  the  number  of 
flight  test  hours  required  for  establishing  the  flutter 
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envelope  of  an  aircraft.  However,  two  fundamental 
issues  must  be  addressed  before  a  method  for  the 
continuous  parametric  identification  of  an  accelerating 
aircraft  can  be  developed. 

The  first  issue  deals  with  determining  how  the 
aeroelastic  properties  of  an  aircraft  can  be  affected 
by  a  constant  acceleration  in  a  level  flight  or  during 
a  maneuver.  In  particular,  it  is  important  to  figure, 
out  whether  it  is  possible  to  relate  in  a  simple  manner 
the  aeroelastic  parameters  measured  in  an  accelerated 
flight  to  those  measured  in  stabilized  flight  conditions. 
To  the  best  of  our  knowledge,  this  issue  has  not  been 
addressed  in  depth  in  the  literature  [1]. 

The  second  issue  is  related  to  the  fact  that  most 
if  not  all  identification  methods  used  in  practice 
implicitly  assume  that  the  given  aeroelastic  system 
is  linear  and  time-invariant.  Whether  these  methods 
can  still  be  used  to  analyze  accelerated  flight  data,  or 
whether  new  methods  are  required  for  this  purpose 
remains  an  open  question. 

In  this  paper,  our  main  objective  is  to  address  some 
preliminary  aspects  of  the  two  issues  raised  above  by 
performing  both  simple  analytical  investigations,  and 
complex  numerical  simulations.  For  this  purpose,  we 
first  derive  the  equations  of  motion  oi  a  typical  wing 
section  set  into  an  accelerated  motion,  and  assess  the 
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effect  of  the  acceleration  on  the  stiffness  of  this  system. 
This  allows  us  to  address  the  first  issue  raised  above. 
Next,  we  overview  a  computational  framework  for  the 
simulation  of  the  transient  aeroelastic  response  of  an 
accelerating  structure,  and  describe  a  “windowing  ’’ap¬ 
proach  for  identifying  the  parameters  of  a  time- varying 
system.  In  order  to  address  the  second  issue  raised  in 
this  paper,  we  perform  several  numerical  simulations 
for  an  F-16  aeroelastic  configuration.  More  specifi¬ 
cally,  we  determine  a  first  damping  vs.  Mach  number 
curve  by  performing  a  series  of  aeroelastic  computa¬ 
tions  corresponding  to  a  sequence  of  stabilized  flight 
conditions. '  Then,  we  perform  a  single  simulation  of 
the  aeroelastic  response  of  the  F-16  airframe  to  per¬ 
turbations  generated  during  an  accelerated  flight.  We 
process  the  continuous  signal  generated  by  this  simu¬ 
lation  by  our  window-based  identification  method,  and 
extract  a  second  damping  vs.  Mach  number  curve  for 
this  aircraft.  We  compare  both  simulated  curves,  con¬ 
trast  them  with  flight  test  data,  and  offer  preliminary 
conclusions  regarding  the  feasibility  of  extracting  the 
flutter  envelope  of  an  aircraft  from  accelerated  flight 
test  data. 

2  Aeroelastic  response  of  a  typical 
wing  section  immersed  in  an 
accelerating  flow 

2.1  Discrete  structural  equations  of  motion 

In  the  absence  of  any  acceleration,  the  equations  of 
motion  of  the  typical  wing  section  shown  in  Fig.  1  can 
be  written  as  [1] 

mh  +  Se8  +  Khh  =  Qh  (1) 

S9h  +  Io8  +  Ke6  =  Qe  (2) 

/ 

whe^e  »ft  is  the  total  mass  per  unit  span  of  the  wing, 
h  t^iid  9  are  the  bending  and  torsional  degrees  of 
freedom  (d.o.f.)  of  the  wing  section,  So  =  m{xQ  —  %c) 
and  Iq  are  its  static  and  polar  moments  of  inertia 
about  the  elastic  axis,  G  and  C  designate  its  center  of 
gravity  and  elastic  center,  Kh  and  K$  are  the  bending 


Fig.  1  A  typical  wing  section:  elastic  center  (C), 
center  of  gravity  (G),  fuselage  (F),  angle  of  attack 
(a),  free-stream  velocity  (14o),  a  typical  point  in 
the  flow  domain  (N) 


and  torsional  stiffness  coefficients,  and  Qh  and  Qo  are 
the  aerodynamic  resultant  force  and  moment  acting 
on  this  section. 

Let  F  denote  a  fixed  point  on  the  fuselage  of  the 
aircraft  whose  typical  wing  section  is  shown  in  Fig.  1. 
If  the  aircraft  accelerates  and  9  remains  small  during 
the  flight,  the  dynamic  equations  of  equilibrium  of  the 
typical  wing  section  become 

mh  +  Sq9  +  Khh  =  Q*h  (3) 

Soh  + 1&9  +  K#9  =  Q*0  (4) 

where 


KMg  =  Ke-  mjFl  (x°G  -  x°c) 

-  miFv  {y°a  -  y°c) 

(5) 

Ql  =  Qh~  mjFv 

(6) 

Qe  =  Qe  +  rn'yp,  (; y%  -  y%) 

~  miFy  {x%  -x%) 

(7) 

and  7fv  are  the  horizontal  and  vertical  compo¬ 
nents  of  the  acceleration  field,  respectively,  and  the 
superscript  0  is  used  to  designate  the  value  at  time 
t  =  0.  Note  that  a  positive  7 fz  corresponds  to  a 
deceleration  of  the  wing  (see  Fig.  1). 

From  Eq.  (6)  above,  it  follows  that  accelerating  the 
wing  section  introduces  not  only  additional  inertia 
forces,  but  also  modifies  the  torsional  stiffness  of  the 
wing  section.  However,  the  reader  can  check  that 
for  level  accelerated  flight  (i.e.  =  0)  and  for 

most  if  not  all  aircraft  configurations,  the  quantity 
7717 f*(xg  "  xc)  IS  Seneral  negligible  compared  to 
Kq .  For  this  reason,  we  reasonably  argue  that  for 
all  practical  purposes,  K#  &  K#.  Since  all  other 
coefficients  of  the  left  hand-side  of  Eqs.  (3,4)  are  not 
affected  by  the  aircraft  acceleration,  we  conclude  that, 
at  least  for  an  accelerated  level  flight,  the  aeroelastic 
parameters  (i.e.  frequency  and  damping  values)  of 
a  typical  wing  section  are  the  same  as  for  stabilized 
flight  conditions. 

In  all  numerical  simulations  that  we  discuss  later  in 
this  paper,  we  time-integrate  Eqs.  (3,4)  by  the  mid¬ 
point  rule  (or  Newmark  method  with  p  =  1/4  and 
7=1/2). 

2.2  Discrete  fluid  equations  of  motion 

The  complete  description  of  the  typical  wing  sec¬ 
tion  introduced  above  requires  the  evaluation  of  the 
aerodynamic  resulting  force  and  moment,  Qh  and  Qo- 
Rather  than  adopting  a  linear  flow  theory  whose  limi¬ 
tations  are  well-known,  we  choose  for  this  purpose  to 
model  the  fluid  by  the  Euler  equations.  More  specif¬ 
ically,  we  rely  on  the  Arbitrary  Lagrangian  Eulerian 
(ALE)  form  [2, 3]  of  the  Euler  equations  in  order  to 
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handle  the  vibrations  of  the  airfoil  and  its  acceleration. 

Let  fl(t)  C  TZ2  be  the  flow  domain  surrounding 
the  airfoil  of  the  typical  wing  section,  and  T(t)  be  its 
moving  boundary.  We  introduce  a  mapping  function 
between  Cl(t)  where  time  is  denoted  by  t  and  a  grid 
point’s  coordinates  by  x,  and  a  reference  configuration 
fi( 0)  where  time  is  denoted  by  r  and  a  grid  point’s  co¬ 
ordinates  by  £,  as  follows 


where  A”  =  ^(Xn),  X2+1  and  X"*1  are  the  follow¬ 
ing  averaged  position  and  velocity  of  the  fluid  moving 
mesh 


X*}  =  a„+1(™^)-^(^-s 

yn+i  A'~+ ' - y  ■ 

Aava  ”  an+lV  Atn  '  <  V  Atn~l 


Kavg 


) 

)  (13) 


an_i,  an,  and  an+i  axe  time-dependent  because  a 
variable  time-step  At"  is  employed  and  are  given  by 


x  =  x(£,r);  t  =  r  (8) 

The  ALE  conservative  form  of  the  Navier-Stokes 
equations  describing  viscous  fio^vs  on  dynamic  meshes 
can  be  written  as  [3-5] 

+  JVx.Fc(W,x)=0  (9) 


c 


<*n 


At" 
At"-1  ’ 


0(1-1  = 


-1  -  C)  Qn+1  — 


c2 

i  +  c 

l  +  2< 

i  +  c 


(14) 


It  remains  to  specify  how  Xn,  the  vector  position  of 
the  fluid  grid  points,  is  updated  at  each  time-station 
t". 


FC{W  ,x)  =  F(W)-xW  (10) 

where  a  dot  superscript  designates  a  time  derivative, 
dx 

J  =  det(dx/d£),  x  =  —  \$  is  the  velocity  of  the 
dynamic  fluid  mesh,  VV  is  the  fluid  state  vector 
(conservative  variables),  and  Tc  denotes  the  ALE 
convective  fluxes. 

We  semi-discretize  Eq.  (9)  on  a  triangulation  from 
which  we  derive  a  dual  mesh  defined  by  time- 
dependent  control  volumes  or  cells  Ci(£).  We  resolve 
the  ALE  convective  fluxes  by  a  suitable  Riemann 
solver  [6-9].  The  resulting  semi-discrete  equation  of 
equilibrium  of  the  fluid  is 

jt(AiWi)  +  Fi(W,X,X)=  0  (11) 

where  -  dxSl,  Wi  de^es  the  average  value 

of  W  over  the  cell  Ci(t),  Fi  denotes  the  semi-discrete 
ALE  convective  flux,  W  is  the  vector  formed  by 
the  collection  of  Wi,  and  X  is  the  vector  of  time- 
dependent  grid  point  positions.  Various  expressions 
of  the  flux  approximation  F{(W,  X,  X)  can  be  found 
in  [6-9]. 


For  this  purpose,  we  first  note  that  for  a  typical  wing 
section  problem,  the  simplest  strategy  for  updating  the 
position  of  the  fluid  mesh  is  to  move  it  rigidly  with 
the  airfoil.  Then,  we  also  note  that  the  motion  of  the 
airfoil  is  completely  determined  by  the  motion  of  the 
fuselage  point  F ,  and  the  vibrations  of  the  bending 
and  torsional  d.o.f.  h  and  0 .  It  follows  that  at  each 
time-station  tn ,  the  position  of  any  fluid  grid  point  N 
(see  Fig.  1)  is  given  by 

=  x p  4-  ( x°N  -  xqf)  cos  0n 
-(y0N-y°F)s  inS" 

vn  =  vf  +  +  ( y%  ~  vf ) cos 

+  (x°N  -  x°F )  sin  9n  (15) 


where  hn  and  9n  are  determined  by  solving  the  coupled 
fluid/structure  equations  of  motion  (3,4,11),  and  the 
instantaneous  position  (z£,y£)  of  the  fuselage  point 
F  is  deduced  from  the  specified  acceleration  of  the 
aircraft.  For  example,  for  a  constant  acceleration,  the 
instantaneous  position  of  the  point  F  is  given  by 


Xp 

—  i  TX 

2 If 


2  * 
2 


+ 

+  v°f 


(16) 


We  time-integrate  the  semi-discrete  equations 
of  flow  motion  (11)  on  dynamic  meshes  using  the 
generalized  second-order  implicit  backward  difference 
scheme  developed  in  [10, 11].  This  scheme  satisfies 
the  second-order  discrete  geometric  conservation  law, 
and  retains  second-order  time-accuracy  on  moving 
grids  [10,11].  It  can  be  written  in  compact  form  as 
follows 

an+1A”+1W?+1  +anA?W?  +  an-1A?-1Wr1 
+  A  tnFi(Wn+1,X2+l,X£g1)  =  0  (12) 


Hence,  the  acceleration  of  the  typical  wing  section  is 
transmitted  to  the  fluid  by  Eqs.  (16),  and  accounted 
for  by  the  additional  convection  term  xW  that 
characterizes  the  ALE  form  (9)  of  the  flow  equations. 

Remark.  At  the  first  glance,  the  reader  may  think 
that  Eqs.  (16)  are  missing  the  terms  VoqJ71  and  V^f1, 
where  Vqo  is  the  ffee-stream  velocity.  However,  these 
terms  are  not  missing.  They  are  automatically  taken 
into  account  by  the  initial  conditions  of  the  CFD  simu¬ 
lation  through  the  specified  ffee-stream  Mach  number 

Moo. 
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2.3  Solution  of  the  coupled  fluid/structure 
equations  of  motion 

Throughout  this  paper,  we  solve  the  coupled 
fluid/ structure  discrete  equations  of  motion  by  the 
second-order  time-accurate  staggered  and  subiteration 
free  algorithm  described  in  [12].  As  stated  earlier,  we 
equip  this  staggered  algorithm  with  the  midpoint  rule 
as  a  structural  time-integrator,  and  the  generalized 
second-order  implicit  backward  difference  scheme  (12) 
developed  in  [10, 11]  as  a  flow  time-integrator. 

3  System  identification  using  the  ERA 

3.1  Identification  on  a  window-by- window  basis 

Depending  on  severed  factors  among  which  the  flow 
regime,  an  aeroelastic  system  can  behave  linearly 
or  non-linearly.  This  raises  a  first  albeit  minor 
concern  as  to  the  applicability  of  several  popular 
signal  processing  techniques  to  the  identification  of 
the  aereolastic  parameters  of  an  aicraft,  or  a  typical 
wing  section,  particularly  in  the  transonic  regime. 
Furthermore,  an  accelerated  aeroelastic  system  is  also 
a  time- varying  system.  This  is  essentially  because  the 
mass,  damping,  and  stiffness  properties  of  the  “wet” 
structure  are  influenced  by  the  free-stream  velocity 
Voc  of  this  structure,  and  varies  with  time  during 
an  accelerated  flight.  This  raises  a  second  concern  as 
to  the  applicability  of  these  signal  processing  methods 
to  the  continuous  parametric  identification  of  an 
accelerating  aircraft.  Both  concerns  can  be  addressed 
by  the  windowing  approach  described  below,  which 
is  in  principle  applicable  to  several  identification 
methods. 

First,  we  note  that 

•  some  methods  are  capable  of  identifying  the  fre¬ 
quency  and  damping  coefficient  of  the  lowest 
mode  of  a  structure  using  as  few  as  2  cycles  of 
the  response  of  this  structure  to  an  input  pertur¬ 
bation. 

•  for  a  mode  at  10  Hz,  2  cycles  correspond  to  a  time 
interval  of  0.2  second. 

•  at  an  altitude  of  3,000  meters,  a  level  flight  ac¬ 
celeration  of  0.05  Mach/second  correponds  to  a 
horizontal  acceleration  of  1.6  gs.  Such  an  acceler¬ 
ation  is  beyond  the  reach  of  most  if  not  all  aircraft. 

•  during  a  time-interval  of  0.2  second,  the  speed 
of  an  aircraft  accelerating  at  0.05  Mach/second 
varies  by  0.01  Mach. 

•  a  1%  variation  of  a  Mach  number  M &  has  a 
negligible  effect  on  the  frequency  and  damping  co¬ 
efficient  of  a  wet  structure  cruising  at  M^. 

From  the  above  observations,  we  conclude  that 
within  a  time-window  of  the  order  of  0.2  second,  an 


aeroelastic  system  with  a  first  wet  mode  at  10  Hz 
can  be  considered  for  all  practical  purposes  as  a  time- 
invariant  system.  Furthermore,  within  a  window  of 
that  size,  the  aereolastic  response  of  such  a  system  can 
be  assumed  to  be  linear,  as  long  as  the  structure  itself 
behaves  linearly,  which  is  usually  the  case  for  an  air¬ 
craft  excited  by  an  input  to  a  control  surface.  Hence, 
we  also  conclude  that  at  least  in  principle,  it  should  be 
possible  to  expand  the  flutter  envelope  of  an  aircraft 
using  its  continuous  vibration  response  to  input  per¬ 
turbations  during  and  accelerated  flight  by  applying 
the  following  simple  procedure 

1.  locate  the  time  instances  at  which  the  Mach  num¬ 
bers  of  interest  are  reached  by  the  aircraft. 

2.  for  each  Mach  number  of  interest,  define  a  time- 
window  of  width  equal  approximately  to  0.2  sec¬ 
ond  or  two  cycles  of  the  expected  lowest  frequency. 

3.  within  each  time- window,  apply  a  suitable  signal 
processing  method  to  identify  the  frequency  and 
damping  coefficient  of  the  wet  structure  and  asso¬ 
ciate  these  with  the  Mach  number  for  which  this 
window  is  defined. 

The  window-based  identification  approach  sum¬ 
marized  above  can  be  performed  either  real-time, 
or  off-line  as  a  signal  post-processing.  In  this  work, 
we  choose  to  base  it  on  the  Eigenvalue  Realization 
Algorithm  (ERA)  [13]. 

The  ERA  is  an  identification  method  for  linear  and 
time-invariant  systems.  It  implicitly  assumes  that  the 
dynamic  response  of  the  given  system  is  sampled  at  a 
constant  rate.  It  can  handle  multi-d.o.f.  systems,  and 
is  less  sensitive  to  noise  than  the  logarithmic  decay 
method.  In  order  to  keep  this  paper  self-contained, 
we  overview  next  the  ERA  and  its  fast  implementa¬ 
tion  FastERA  [14, 15]  which  is  suitable  for  real-time 
processing. 

3.2  Overview  of  the  ERA 

In  a  linear  context,  damped  structural  vibrations 
are  governed  by  the  following  equations  of  dynamic 
equilibrium 

Mq  +  Dq  +  Kq  =  f(t)  (17) 

where  M,D,K  are  respectively  the  mass,  damping 
and  stiffness  matrices,  and  where  f  denotes  the  vec¬ 
tor  of  external  loads.  The  above  equations  can  be 
re-written  in  state-space  form  as  follows 

r 

X  =  /Cx  +  £f(f) 

J  (18) 

q  =  Cx 
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where 


Hence,  q*  is  the  classical  solution  of  the  state-space 
realization  problem  for  the  following  linear  sampled 
data  system 


(19) 

where 


f 

Xfc+i  =  Ax*  +  Bu* 

< 

qj t  =  Cx*  +  Gu* 
v. 


(25) 


Assuming  that  initially  the  structural  system  is  at  rest, 
and  that  it  is  excited  at  t  =  0  by  an  impulse  load 
f (t)  =  uo<f(f)  where  5  denotes  the  Dirac  delta  function, 
the  solution  of  Eq.  (18)  is 

q(t)  =  CeKtBuo  (20) 

Note  that  Eq.  (20)  is  also  the  solution  of  the  prob¬ 
lem  (18)  for  f(£)  =  0  and  an  initial  velocity  equal  to 
M-1uo. 

For  an  arbitrary  excitation,  f(i)  can  be  represented 
by  a  series  of  impulses  at  discrete  time-stations  U  — 

oo 

that  is,  f(f)  =  u*£(t-- 1»),  in  which  case  the  solution 

of  Eq.  (18)  at  time  tk  is  given  by 


B  =  A  B  C=C  (26) 

If  ridof  denotes  the  total  number  of  d.o.f.  of  the  com¬ 
putational  model  described  by  Eq.  (17),  but  only 
n£e;  <  ridof  d.o.f.  are  measured,  then  the  matrix 
C  can  be  written  as 

c  =  [  L  (27) 

where  L  is  an  x  ridof  Boolean  matrix.  Hence,  in 
general,  the  number  of  state-space  variables  xk  that 
are  considered  is  nx  =  2rid0/,  and  the  matrix  A  defined 
in  Eq.  (22)  is  an  nx  x  nx  matrix. 

In  general,  the  discrete  convolution  sum  (23)  is  ex¬ 
pressed  as 

k 

qJfe  =  ^Xfc_iUi  (28) 

»=o 


k 

q*  =  q  (tk)  = 

t=0 

=  YlCe{k~i)AtKB\ii  (21) 

i— 0 

where  a  constant  sampling  rate  1/A t  is  assumed  so 
that  tk  -  t{  =  (k  -  i)At. 


where  the  matrices 


t 

G  ,  m  =  0 

< 

CAm-1B  ,  m  >  0 


(29) 


are  known  as  the  Markov  paurameters.  Hence,  can 
be  directly  related  to  the  load  components  u*  via  the 
the  Markov  parameters 


Let 

A  =  eAtlc  (22) 

Using  the  above  definition  of  the  matrix  A,  Eq.  (21) 
can  be  re-written  as  follows 

q*  =  ^CA^Bm 

i=0 

fc-l 

=  ^  CAfc-’-1  (AH)u,-  +  CBuk 

t=0 

k-1 

=  ^(CAfc-’-1Bui)  +  Gu*  (23) 
1=0 

where 

G  =  CB  (24) 


Let  za  and  \s  =  as  ±  ius  denote  the  complex  eigen- 
modes  and  eigenfrequencies  of  the  structure  defined 
by 

X23Mz3  +  A5Czs  .+  Kzs  =  0  (30) 

From  Eqs.  (19),  it  follows  that  the  \$  aire  also  the 
eigenvalues  of  £,  and  the  eigenvectors  of  this  matrix 
are 


The  ERA  exploits  the  results  summarized  above  as 
follows.  First,  it  constructs  the  Markov  parameters 
Mk  using  the  input  and  output  data,  namely,  the 
impulse  loads  and  measured  displacements.  Then,  it 
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extracts  the  matrix  A  from  the  Markov  parameters 
Mk  using,  for  example,  the  Fast  ERA  algorithm  de¬ 
scribed  in  Section  3.3.  Finally,  the  ERA  computes  the 
eigenvalues  o3  and  eigenvectors  aa  of  A  which  satisfy 

Aas  =  a8  a5  (32) 

From  Eq.  (22),  it  follows  that  the  sought-after  complex 
modes  zs,  damping  ratios  and  frequencies  ojs  of  the 
given  structural  system  are  given  by 


Z  S 

=  Ca, 

(33) 

(34) 

Vs 

=  ^Im{ln(CTs)} 

(35) 

3.3  The  FastERA  algorithm  [14] 

In  this  work,  we  use  the  FastERA  implementation 
of  the  ERA  to  perform  all  ^identifications.  For 
this  reason,  we  review  FastERA  in  this  section  and 
introduce  some  key  variables  that  appear  in  the 
remainder  of  this  paper. 


The  FastERA  method  is  based  on  the  analysis  of 
the  following  Hankel  matrix  that  is  defined  for  a  data 


set  (qfc,Ufc)  sampled 

at  N  points  in  time 

All 

M.2 

Md 

8 

•Q 

PL 

II 

M.2 

M3  ■■ 

Md+i 

(36) 

t 

* 

_ i 

Mg+2  ■  ■ 

■  Mq+d-l 

and  where  q  and  d  are  such  that  N  =  q  +  d.  From  the 
definition  (29)  of  the  Markov  parameters,  it  follows 
that  H qd  can  be  decomposed  as  follows 

Hqd  =  V,Wd  (37) 

where 


V,= 


C 

CA 

CA9“: 


,  Wd  = 


B  AB  A^B 


(38) 


The  qn™oy  x  qn™*f  matrix  Vq  determines  the  observ¬ 
ability  of  the  system  whereas  Wd  is  related  to  its 
controllability. 


The  FastERA  method  starts  with  the  construction 
of  the  qn^0ef  x  qnfiff  “square  data  matrix” 

J,  =  H,dH^  =  V,WdWjV^  (39) 

6  OF 


For  a  given  set  of  inputs  u*  and  outputs  q*,  the  choice 
of  state-space  variables  x  is  not  unique.  Indeed,  the 
input /output  relation  is  not  affected  by  the  change 
of  variable  x  =  Tx  which  transforms  the  realization 
system  (25)  into 

,  **«  “  T-‘AT*,+T-‘Bu.  (4o) 

q*  =  CTx*  +  Gu* 

The  controllability  factors  associated  with  x  are 

WdWTd  =  T~1WdWdT~T  (41) 

In  particular,  choosing  T  =  (W^Wj)1^2  leads  to 

-  T 

WdWd  =  I.  Hence,  without  any  loss  of  generality, 
the  FastERA  algorithm  assumes  that  the  state-space 
variables  are  chosen  so  that  Eq.  (39)  simplifies  to 

J,  =  V9V^  (42) 

From  Eq.  (38),  it  follows  that  the  first  (q  - 
rows  of  Vq  are  given  by 


V<i)  = 


c 

CA 


CA9 


-2 


(43) 


and  the  last  ( q  -  l)n^ys  rows  of  that  matrix  are  given 
by 


VW  = 


CA 

CA9"1 


=  Vj‘)A 


In  general,  the  number  of  d.o.f.  ndo/  associated  with 
a  system  to  be  identified  is  unknown  a  priori.  Indeed, 
ridof  is  infinite  for  any  continuous  system.  In  prac¬ 
tice,  the  target  values  for  n™*fS  and  the  corresponding 
number  of  state-space  variables  nx  =  2 ndo/  are  dic¬ 
tated  by  the  complexity  of  the  model  to  be  realized. 
Once  these  values  are  set  by  the  user,  the  FastERA 
algorithm  computes  the  nx  largest  eigenvalues  kj  of 
J q  and  their  corresponding  eigenvectors  pj  in  order  to 
build  the  square  root  factorization 

J,  =  Pdiag(«j)PT  (45) 

=  (Pdiag(«;j/2))  (diag(/c]/2)PT) 

=  V,V[  (46) 
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where  P  =  [p:  ...  ...]  and  V,  =  Pdiag(K;V2).  From 

Eqs.  (43,44)  applied  to  V,  it  follows  that 

A=(vi1))+v;2)  (47) 

where  is  the  pseudo-inverse  of 

Hence,  the  extraction  of  the  nx  largest  eigenpairs 
of  the  x  qn™*j  data  matrix  Jq  is  the  most 

computationally  significant  step  of  the  FastERA 
algorithm.  Keeping  q  relatively  small,  say  q  w  d/5, 
allows  FastERA  to  operate  in  real  time.  For  further 
details  on  this  identification  method,  we  refer  the 
reader  to  [14]. 

Next,  we  validate  the  window-by-window  applica¬ 
tion  of  the  ERA  for  the  indentification  of  the  parame¬ 
ters  of  an  accelerating  aircraft  by  a  series  of  numerical 
simulations  designed  for  an  F-16  fighter  configuration 
for  which  flight  test  data  is  available. 

4  Simulation  of  the  continous 
parametric  identification  of  an  F-16 
configuration 

Here,  we  focus  on  a  Block  40  version  of  the  F-16 
fighter.  Because  of  CPU  time  limitations,  we  consider 
only  two-dimensional  numerical  simulations.  There¬ 
fore,  we  report  first  on  the  design  of  a  typical  wing 
section  model  for  this  aircraft. 

4.1  A  typical  wing  section  model  for  the  F-16 

We  start  from  a  detailed  finite  element  structural 
model  of  a  “clean”  right  wing  of  the  F-16  Block  40 
equipped  with  a  missile  launching  system  at  its  tip 
(Fig  2).  Our  objective  is  to  construct  a  two-d.o.f. 
wing  section  model  that  is  equivalent  to  the  three- 
dimensional  wing  in  the  following  sense 

1.  it  reproduces  the  first  bending  and  torsion  fre¬ 
quencies  which  are  predicted  by  the  three- 
dimensional  finite  element  model  of  the  wing  to 
be  4.76  Hz  and  7.43  Hz,  respectively. 

2.  it  reproduces  the  same  bending  and  torsion  modal 
masses. 

3.  it  reproduces  the  same  vertical  displacement  at 
the  leading  edge  of  the  cross  section  located  at 
68%  of  the  distance  from  the  root  to  the  tip  of 
the  three-dimensional  wing,  for  both  bending  and 
torsion  mode  shapes,  when  these  are  normalized 
to  a  unit  rotation  of  the  cross  section. 

When  the  wing  has  a  large  aspect  ratio  and  a  small 
sweep  angle,  it  is  commonly  suggested  to  choose  the 
mechanical  properties  of  the  typical  wing  section  as 
to  match  those  of  the  cross  section  located  at  70%  to 


75%  of  the  distance  from  the  root  to  the  tip  of  the 
three-dimensional  wing  (for  example,  see  [1]  and  ref¬ 
erences  6-1  and  6-2  therein).  Nevertheless,  we  choose 
in  criterion  3  stated  above  the  cross-section  located  at 
68%  of  the  distance  from  the  root  to  the  tip  (see  Fig. 
2)  because  the  F-16  wing  is  strongly  tapered  and  is 
rather  soft  towards  its  tip. 

Besides  the  shape  of  the  airfoil,  six  parameters  de¬ 
fine  the  sought-after  typical  wing  section,  These  are 
denoted  collectively  by  Ptws ,  where  the  subscript  tws 
stands  for  typical  wing  section,  and  listed  below 

Ptws  =  {mJo,xG,  {xq  -  xc),KhiKo} 


The  three  criteria  stated  above  for  establishing  the 
equivalence  between  the  typical  wing  section  and  the 
three-dimensional  wing  can  be  formulated  as  follows 


^tW3,l{Ptws) 

2tt  x  4.76  rad/s 

^•tws^^P ttys) 

27 r  x  7.43  rad/s 

{Ptws) 

1.375  106  Kg.m2 

l*twSy2{Ptws) 

2.523  Kg.m2 

^tws^li^tws) 

25.61  m 

titeWSy2(Ptws) 

1.017  m 

where  ft  denotes  a  frequency,  the  subscripts  1  and  2 
refer  to  the  bending  and  torsion  modes,  respectively, 
/i  denotes  a  modal  mass,  the  superscript  le  designates 
the  leading  edge,  hle  =  h  4-  ( xle  -  xc){6  —  1)  is 
the  vertical  displacement  at  the  leading  edge  of  the 
typical  wing  section,  and  all  components  of  the  right 
hand-side  of  Eq.  (48)  are  obtained  from  the  detailed 

2 


1:  cross-section  chosen  for  the  structural  properties 
2:  cross-section  chosen  for  the  aerodynamic  properties 

Fig.  2  Three-dimensional  detailed  finite  element 
model  of  an  F-16  wing 
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finite  element  model  of  the  three-dimensional  wing. 
Note  that  the  fact  that  hlteW8jl  »  hlteWSy2  suggests  that 
despite  the  bending/torsion  coupling,  the  first  mode 
of  the  wing  is  dominated  by  bending,  while  the  second 
one  is  dominated  by  torsion. 

The  constraints  (48)  lead  to  a  nonlinear  system  of 
equations  with  six  unknowns,  which  we  solve  by  the 
Nelder-Mead  simplex  method  (function  fminf  in  Mat- 


lab)  to  obtain 

m 

=  2.05  103  Kg 

Ie 

=  2.53  10JKg.m2 

XQ 

=  1.128  m 

(upstream  of  section  at  68%) 

(xq  -  Xc ) 

=  0.0642  m  (G  behind  C) 

Kh 

=  2.046  106  N/m 

Ke 

=  5.468  106  N.m 

Using  the  above  numbers,  we  find  that  the  modified 
torsional  stiffness  (see  Eq.  6)  for  an  acceleration  as 
high  as  0.05  Mach/second  is  K #  =  5.471  106  N.m. 
Hence,  as  expected,  K$  «  Kq,  which  supports  the 
conclusion  made  in  Section  2.1  as  to  the  little  effect  of 
level  flight  acceleration  on  the  aeroelastic  parameters 
of  a  system. 

Finally,  we  assume  modal  damping  ratios  of  1.0% 
and  1.5%  for  the  first  and  second  vibration  modes  of 
the  typical  wing  section,  respectively. 


4.2  The  airfoil  and  the  flow  computational 
domain 

We  choose  the  airfoil  of  .he  typical  wing  section  as 
that  of  the  cross  section  located  at  45%  of  the  distance 
between  the  root  and  tip  of  the  three-dimensional  wing 
(see  Fig.  2),  because 

•  the  chord  of  the  airfoil  of  the  typical  wing  section 
must  be  close  to  the  ratio  of  the  wetted  area  and 
the  wing  span. 

•  because  of  tapering,  most  of  the  lift  is  generated 
by  the  section  of  the  wing  that  is  close  to  the  root, 
which  means  that  the  aerodynamic  center  of  the 
wing  is  within  that  area. 

We  discretize  the  flow  domain  by  18,000  vertices, 
and  ensure  a  sufficient  resolution  for  shock  capturing 
in  the  region  close  to  the  sharp  leading  edge  (Fig.  3). 
Because  the  purpose  of  the  typical  section  is  to  rep¬ 
resent  the  entire  wing,  we  multiply  each  areodynamic 
force  obtained  from  a  flow  computation  on  this  two- 
dimensional  mesh  by  the  span  of  the  wing. 


Fig.  3  Discretization  of  the  flow  computational 
domain  (F-16  airfoil) 

4.3  Aeroelastic  numerical  simulations 

We  fix  the  altitude  to  3,000  m,  set  the  angle  of 
attack  to  0  degree,  and  perform  first  a  series  of  aeroe- 
latic  simulations  for  a  sequence  of  stabilized  flight 
conditions  at  the  following  Mach  numbers:  0.8,  0.85, 
0.875,  0.9,  0.925,  0.95,  1.0,  1.1,  1.2,  1.3  and  1.4.  For 
each  Mach  number,  we  start  the  numerical  simulation 
from  a  steady-state  flow  and  the  following  initial 
conditions  for  the  typical  wing  section:  h°  =  0.01 
m/s,  and  9°  =  0.2  rad/s.  We  report  in  Fig.  4  and  Fig. 
5  the  predicted  transient  responses  of  the  structure 
at  Mqo  =  0.8  and  Moo  =  1.0,  respectively.  We  also 
include  at  the  top  of  both  figures  the  iso-Mach  contour 
plots  at  t  =  0.  These  figures  show  that  at  =  0.8, 
the  aeroelastic  vibrations  are  rapidly  damped  out. 
They  also  show  that  at  M oo  =  0.8,  both  modes  of 
the  typical  wing  section  contribute  initially  to  the 
bending  d.o.f.  h,  but  only  one  mode  of  the  typical 
wing  section  contributes  to  the  history  of  the  torsional 
d.o.f.  9.  We  conclude  that  it  is  the  bending  mode  of 
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the  typical  wing  section  that  is  rapidly  damped  out  at 
Moo  =  0.8.  The  history  of  9(t)  graphically  depicted 
in  Fig.  5  reveals  that  flutter  or  limit  cycle  oscillations 
are  initiated  at  Moo  =  1*0. 

Next,  we  simulate  the  aeroelastic  response  of  the 
typical  wing  section  during  a  flight  accelerated  at  the 
rate  of  0.05  Mach/second.  We  remind  the  reader  that 
such  an  acceleration  is  even  higher  than  what  an  F-16 
can  achieve  in  a  level  flight.  We  initiate  this  simula¬ 
tion  at  Moo  =  0.75  and  excite  the  structure  with  the 
same  initial  conditions  as  previously:  h°  =  0.01  m/s, 
and  6 0  =  0.2  rad/s.  We  report  in  Fig.  6  the  predicted 
response  of  the  structure.  The  reader  can  observe  that 
both  modes  of  the  typical  wing  section  contribute  to 
the  early  response  of  the  h  d.o.f.  The  significant  de¬ 
crease  of  the  mean  value  of  h  that  is  noted  between 
Mach  0.8  and  Mach  0.95  suggests  an  important  drop  in 
the  lift  in  this  Mach  region,  which,  given  that  the  F16 
airfoil  is  unsymmetric,  is  indicative  of  the  appearance 
of  a  shock  and  reaching  the  Mach  divergence  speed. 
Most  importantly,  Fig.  6  shows  that  between  Mach 
0.85  and  Mach  0.95,  the  vibrations  of  h  and  6  become 
too  small  to  allow  a  parametric  identification.  Hence, 
we  conclude  that  a  continuous  parametric  identifica¬ 
tion  of  an  aicraft  that  accelerates  across  the  subsonic, 
transonic,  and  supersonic  regimes  requires  systematic 
re-excitations. 

In  order  to  illustrate  the  effect  of  re-excitation,  we 
simulate  a  second  accelerated  flight  where  the  initial 

JpnS  ,  ,  . *  ■  *  <  w, ,  U  f.M'i't  *'i  f.fygl 

( ,  '  ,  ’  '  ‘  ‘  ' ,  ,1  , , ,  *  *  ,  tjjj 

aE  ‘  , ,  ,  ’  i  ,  ' , ' "  ‘  '  'I  pi  i.i***» 

f  Jp  *  >  *  •  •  •  ,  *  1  .  *  -  *  •’  ’  _  •  Ji 


Fig.  4  Transient  aeroelastic  response  for  stabilized 
flight  conditions  at  Mach  =  0.80 


Fig.  5  Transient  aeroelastic  response  for  stabilized 
flight  conditions  at  Mach  =  1.00 


Fig.  6  Transient  response  during  an  accelerated 
flight:  initial  speed  corresponds  to  Mach  0.75  and 
is  increased  by  0.05  Mach/ second 
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mu 


speed  corresponds  to  Mach  number  0.84,  and  the  ini¬ 
tial  excitation  is  effected  by  the  same  initial  conditions 
as  previously.  We  report  in  Fig.  7  the  computed 
aeroelastic  response.  This  response  is  characterized 
by  larger  amplitudes  of  vibration  between  Mach  0.85 
and  Mach  0.95,  and  confirms  that  flutter  or  limit  cycle 
oscillations  initiate  around  Mach  1.0. 

Finally,  we  report  in  Fig.  8  the  simulated  aeroelastic 
response  of  the  accelerated  typical  wing  section  after 
re-excitation  at  Mach  1.04. 

4.4  Parameter  identification  using  the  ERA: 
stabilized  vs.  accelerated  flight  scenarios 

As  stated  in  Section  3,  the  ERA  assumes  a  constant 
sampling  rate.  Since  our  numerical  simulations  are 
not  performed  with  a  constant  time-step  At,  we 
post-process  our  simulation  results  by  a  quadratic 
interpolation  scheme  in  order  to  generate  signals  with 

h  (m) 

0.02 
0.01 


-0.01 

-0.02 

-0.03 

-0.04 

-0.05 

e  (° 


0.85  0.9  0.95  1  1.05  1.1 

Mach 

Fig.  7  Transient  response  during  an  accelerated 
flight:  initial  speed  corresponds  to  Mach  0.84  and 
is  increased  by  0.05  Mach/second 

o 

h  (m) 

-0.01 

-0.02 

-0.03 

-0.04 

-0.05 

e  (°)0-8 

0.6 

0.4 

0.2 

0 

-0,2 

1.1  1.2  1.3 

Mach 

Fig.  8  Transient  response  during  an  accelerated: 
initial  speed  corresponds  to  Mach  1.04  and  is  in¬ 
creased  by  0.05  Mach/second 


a  constant  sampling  rate. 

The  typical  wing  section  described  in  this  paper  has 
2  d.o.f.  However,  the  aeroelastic  typical  wing  section 
has  more  than  2  d.o.f.,  because  of  the  surrounding 
fluid.  For  this  reason,  in  all  cases  discussed  in  this 
section,  we  set  the  number  of  states  for  the  synthesized 
system  to  nx  =  10,  as  if  the  system  contained  5  d.o.f. 

Our  experience  is  that  the  ERA  requires  two  cy¬ 
cles  of  the  lowest  mode  contributing  to  the  signal,  and 
about  500  sampling  points  per  cycle  in  order  to  iden¬ 
tify  accurately  our  system.  Hence,  in  this  work  we  set 
the  parameters  of  the  ERA  as  follows 


1000  < 

N 

<  6000 

500  < 

Q 

<  1000 

500  < 

d 

<  5000 

2Tl  < 
2ns 

A  ts 

<3ri 

~N 

nx 

- 

10 

where  A ts  is  the  sampling  time,  and  T\  is  the  period  of 
the  lowest  mode.  We  have  verified  a  posteriori  that  the 
ERA  configured  with  the  above  parameters  produced 
excellent  results  for  all  applications  discussed  in  this 
paper.  Nevertheless,  we  note  that  higher  values  of  the 
number  of  samples  N  and  higher  values  of  q  improve 
the  accuracy  of  the  identification,  but  increase  its  cost. 

As  explained  in  Section  3,  in  order  to  extract  the 
aeroelastic  parameters  of  the  typical  wing  section  from 
the  signals  generated  by  the  accelerated  flight  simula¬ 
tions,  we  propose  to  employ  the  ERA  with  windowing. 
The  size  of  each  time- window  is  N  At  s-  Given  that  the 
frequency  of  the  first  wet  mode  of  the  system  can  be 
expected  to  be  close  to  that  of  the  first  dry  mode  of 
the  system  —  that  is.  4.7  Hz  —  it  follows  that  the 
size  of  each  time-window  varies  between  0.4  second 
and  0.6  second.  Hence,  for  an  acceleration  of  0.05 
Mach/second,  the  maximum  variation  of  the  Mach 
number  within  a  time- window  is  about  0.03  Mach. 
The  significance  of  this  variation  depends  on  how  fast 
the  damping  coefficient  of  the  structure  varies  with  the 
Mach  number,  which  depends  on  the  flow  regime.  Ex¬ 
tensive  experiments  have  revealed  that  when  using  the 
ERA  with  windowing,  the  identification  results  are  to 
some  extent  insensitive  to  small  variations  in  the  size 
of  the  time- window  and/or  the  values  of  ns,  q ,  and  d. 

The  frequencies  and  damping  ratios  identified  by 
the  FastERA  using  stabilized  and  accelerated  flight 
simulation  data  are  reported  in  Fig.  9.  The  following 
observations  are  noteworthy 

•  our  typical  wing  section  exhibits  flutter  for  the 
second  (torsional)  mode  in  the  transonic  regime, 
for  0.89  <  Moo  <  1.1. 
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Mode  1  (bending) 

„  frequency  (Hz) _ ^ _ t_ 


fied  because  it  becomes  significantly  damped  out. 
Re-exciting  the  system  every  second  or  so  cures 
this  problem. 

for  the  second  mode,  the  frequencies  and  damping 
ratios  identified  for  the  stabilized  and  accelerated 
flight  scenarios  are  in  good  agreement. 


’  0.0  0.9  1  1.1  1-2  1-3  1.4 


-1 - from  stabilized  flight  conditions 

o  from  accelerated  flight  (0.05  Mach/second) 


Fig.  9  Identified  frequencies  and  damping  ratios 

•  a  sharp  decrease  of  the  damping  ratio  82  of  the 
torsional  mode  occurs  around  Mach  0.9,  and  a 
slow  increase  of  that  damping  ratio  occurs  above 
Mach  0.92.  The  same  trend  is  also  observed  for 
the  damping  ratio  8\  of  the  first  mode;  however, 
8\  remains  positive. 

•  the  frequency  of  the  first  mode  appears  to  be  al¬ 
most  independent  of  the  Mach  number.  On  the 
other  hand,  the  frequency  for  the  torsional  mode, 
which  is  responsible  here  for  flutter,  slightly  in¬ 
creases  with  the  Mach  number. 

•  occasionally,  the  frequencies  and  damping  ratios 
identified  for  the  first  mode  in  simulated  acceler¬ 
ated  flight  are  reported  to  be  different  from  those 
identified  using  stabilized  flight  simulated  data. 
This  is  due  to  the  fact  that  in  our  accelerated 
flight  simulation,  because  the  system  is  excited 
only  at  the  beginning  of  the  flight  segment,  after 
a  certain  amount  of  time,  the  contribution  of  the 
first  mode  to  the  signal  can  no  longer  be  identi¬ 


5  Validation 


Here,  we  compare  our  simulated  flight  results  to 
flight  test  data  provided  by  the  Edwards  Air  Force 
Base.  The  simulation  results  and  test  data  are  for  the 
same  stabilized  flight  conditions,  but  not  the  same  F- 
16  configuration.  The  typical  wing  section  designed  in 
this  paper  is  for  a  clean- wing  configuration  of  the  F-16. 
The  flight  test  data  provided  by  the  Edwards  Air  Force 
Base  is  for  a  configuration  of  the  F-16  that  includes  py¬ 
lons  and  missiles.  Nevertheless,  Fig.  10  shows  that  our 
predicted  wet  frequencies  are  in  good  agreement  with 
those  measured  in  flight  test.  However,  our  predicted 
damping  ratios  do  not  agree  well  with  the  experimen¬ 
tal  data,  even  though  the  trend  of  their  variation  with 
the  Mach  number  is  similar  to  that  of  the  flight  test 
data.  We  can  reasonably  argue  that  a  major  cause 
of  this  discrepancy  is  the  typical  wing  section  model, 
which  is  supposed  to  be  realistic  only  for  fairly  homo¬ 
geneous  wings  with  high  aspect  ratios  and  small  angles 
of  sweep.  Other  possible  causes  are  the  pylons,  mis¬ 
siles,  and  viscous  effects  that  are  not  accounted  for  in 
our  flow  computations.  Perturbing  the  parameters  of 
the  typical  wing  section  can  improve  the  simulation 
results.  For  example,  if  the  position  of  the  elastic  cen¬ 
ter  of  the  typical  wing  section  is  shifted  to  35%  of  the 
chord,  our  numerical  simulations  produce  results  that 
are  in  better  agreement  with  the  flight  test  data  as 
demonstrated  in  Fig.  11. 

Mode  2  (torsional) 


V  flight  test  data 


Fig.  10  Validation 
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Mode  2  (torsional) 


Mach 


— I —  numerical  simulation,  xc  =  6%  of  chord 
— © —  numerical  simulation,  xc  =  35%  of  chord 

V  flight  test  data 

Fig.  11  Sensitivity  of  the  aeroelastic  parameters 
to  the  position  of  the  elastic  center 

6  Conclusions 

The  wet  frequencies  and  damping  ratios  of  an  air¬ 
craft  are  functions  of  the  free-stream  Mach  number. 
Therefore,  during  an  accelerated  flight,  these  parame¬ 
ters  become  time-dependent.  For  this  reason,  theory 
suggests  that  standard  signal  processing  techniques, 
which  are  limited  to  time-invariant  systems,  cannot 
be  applied  to  the  identification  of  an  accelerated  aeroe¬ 
lastic  system.  However,  within  a  time-window  of  0.2 
second,  and  for  a  typical  level  flight  acceleration  of  lg 
or  0.03  Mach/second,  the  Mach  number  varies  by  0.6% 
only.  Hence,  within  a  time-window  of  0.2  second,  the 
aeroelastic  parameters  of  an  aircraft  can  be  assumed  to 
r  main  constant.  Such  a  time- window  corresponds  to  2 
ry:les  of  a  mode  at  10  Hz,  a  frequency  that  is  relevant 
to  the  first  bending  and  torsion  modes  of  most  aircraft. 
Therefore,  any  signal  processing  technique  that  can 
identify  correctly  a  mode  from  2  cycles  of  its  response 
is  a  candidate  technique  for  identifying  the  parameters 
of  an  accelerated  aeroelastic  system  in  a  window-by- 
window  approach.  Furthermore,  the  analytical  study 
of  the  typical  wing  section  shows  that  realistic  level 
flight  accelerations  do  not  affect  the  frequencies  and 
damping  ratios  of  an  aeroelastic  system.  In  other 
words,  the  aeroelastic  parameters  of  an  aircraft  identi¬ 
fied  during  an  accelerated  flight  are  the  same  as  those 
identified  in  stabilized  flight  conditions.  Hence,  we 
conclude  that  there  is  a  hope  that  flutter  testing  could 
be  performed  in  accelerated  flight,  thereby  reducing 
the  cost  and  risk  involved  in  determining  the  flutter 
envelopes  of  fighters.  In  practice,  technical  details 
such  as  re-excitation  procedures,  real-time  identifica¬ 
tion,  and  flutter  early  alert  systems  must  be  addressed 
to  enable  such  a  technology. 
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Abstract 

We  overview  a  three-field  formulation  of  coupled  fluid-structure  inter¬ 
action  problems  where  the  flow  is  modeled  by  the  arbitrary  Lagrangian- 
Eulerian  form  of  either  the  Euler  or  Navier-Stokes  equations,  the  struc¬ 
ture  is  represented  by  a  detailed  finite  element  model,  and  the  fluid  grid 
is  unstructured,  dynamic,  and  constructed  by  a  robust  structure  analogy 
method.  We  discuss  the  latest  advances  in  the  computational  algorithms 
associated  with  this  approach  for  modeling  aeroelastic  problems.  We  ap¬ 
ply  the  three-field  nonlinear  computational  framework  to  the  prediction 
of  the  aeroelastic  frequencies  and  damping  coefficients  of  an  F-16  con¬ 
figuration  in  various  subsonic,  transonic,  and  supersonic  regimes.  We 
consider  for  this  purpose  both  the  popular  two-dimensional  typical  wing 
sr^Uon  model  and  a  detailed  three-dimensional  finite  element  model  of 
the  structure,  and  compare  in  both  cases  the  obtained  numerical  results 
with  flight  test  data.  We  comment  on  the  advantages  and  shortfalls  of 
both  approaches,  and  on  the  feasibility  as  well  as  the  merit  of  the  three- 
field  formulation  of  nonlinear  aeroelasticity  for  the  extraction  of  flutter 
envelopes. 

Key  words.  Aeroelasticity,  F-16  fighter,  fluid-structure  interaction,  flutter,  CFD  on 
moving  grids. 
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1  Introduction 


If  an  elastic  aircraft  immersed  in  an  unsteady  flow  undergoes  a  damped  harmonic 
motion  characterized  by  small  displacement  amplitudes ,  a  wet  (or  aeroelastic) 
circular  frequency  Q,  and  a  wet  positive  or  negative  damping  coefficient  a,  and  if 
the  airflow  surrounding  it  can  be  accurately  predicted  by  an  inviscid  linearized 
theory ,  then  aeroelastic  stability  problems  such  as  flutter  can  be  cast  into  an 
eigenvalue  problem  of  the  form  [1] 

1  2a 

(1) 

where  Zm  €  Rmxm  is  given  by 

zm(k)  =  -  \^Am(k))  (2) 

and  i  is  the  complex  imaginary  number  satisfying  i2  —  —  1.  Here,  Moo, 
and  poo  denote  the  free-stream  Mach  number,  velocity,  and  density  of  the  flow. 
k  =  lj/Voo  is  the  reduced  aeroelastic  circular  frequency  associated  with  u.  Sl2m 
is  a  diagonal  matrix  storing  the  squares  of  the  first  m  fundamental  circular 
frequencies  of  the  dry  aircraft,  and  Im  is  the  identity  matrix  of  dimension 
m.  Typically,  the  aircraft  structure  is  represented  by  its  first  20  to  40  dry  (or 
ground)  modes,  and  therefore  20  <  m  <  40.  Am  denotes  the  projection  of  a 
linear  aerodynamic  operator  A  onto  the  m-dimensional  modal  basis  associated 
with  and  ym  denotes  the  projection  of  the  amplitude  of  the  aircraft’s 
motion  onto  that  basis. 

A  critical  value  of  jfc,  fccr,  is  that  for  which  Zm  has  a  real  eigenvalue  equal 
to  l/£fr2  (1  <  i  <  m).  From  Eq.  (1),  it  follows  that  for  k  =  fccr,  the  damping 
coefficient  a?r  of  the  z-th  aeroelastic  mode  vanishes.  Hence,  given  a  free-stream 
Mach  number  and  a  free-stream  density,  the  critical  value  of  the  reduced  fre¬ 
quency,  if  it  exists,  can  be  found  by  sweeping  on  k  and  solving  the  eigenvalue 
problem  (1)  until  a  real  eigenvalue  is  found.  In  that  event,  the  flutter  mode  of 
the  aircraft  is  the  mode  i  for  which  af  =  0,  the  flutter  speed  of  the  aircraft 
is  V£r  =  Qfr/kcr ,  and  the  flutter  dvnamic  pressure  is  Acr  =  PooV^f  /2.  Such 
a  procedure  for  extracting  the  flutter  speed  of  an  aircraft  is  known  as  the  ak” 
method  [2].  It  is  accurate  when  uhe  assumptions  stated  above  hold,  and  when 
the  structure  is  less  than  10%  damped.  When  the  structure  has  a  higher  per¬ 
centage  of  damping,  an  improved  version  of  this  procedure  known  as  the  “p-k” 
method  [3]  is  preferable  for  finding  the  flutter  dynamic  pressure.  Assuming  that 
the  positive  or  negative  a*  coefficients  are  small,  the  variations  of  the  aeroelas¬ 
tic  frequencies  and  damping  coefficients  with  the  Mach  number,  u^(Moo)  and 
a*  (Moo),  l  <i  <7n,  can  also  be  predicted  by  sweeping  on  M^  and  solving  the 
eigenvalue  problem  (1). 

The  “k” ,  “p-k”  and  other  similar  procedures  for  flutter  analysis  are  based 
on  the  linear  theory  of  aeroelasticity  whose  assumptions  have  been  stated  at 
the  beginning  of  this  introduction.  They  are  fast  and  memory  lean.  For  this 
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reason,  they  are  also  popular  in  the  aerospace  industry.  In  the  subsonic  regime, 
most  if  not  all  of  these  procedures  rely  today  on  the  doublet-lattice  method  [4] 
for  computing  the  linear  aerodynamic  operator  A.  It  is  interesting  to  note  that 
this  method,  which  was  developed  over  thirty  years  ago,  is  still  today  the  most 
used  method  for  predicting  subsonic  unsteady  flows  in  production  environments, 
particularly  for  load  and  flutter  analyses.  In  the  supersonic  regime,  various 
methods  related  to  the  piston  theory  [5]  are  used  for  constructing  A. 

High  performance  military  aircraft  are  usually  flutter  critical  in  the  transonic 
speed  range  at  high  dynamic  pressure.  Unfortunately,  in  that  regime,  the  mixed 
subsonic-supersonic  flow  patterns  and  shock  waves  are  such  th^t  the  linear  flow 
theory  in  general  —  and  therefore  the  doublet-lattice  method  in  particular  — 
axe  not  reliable  for  predicting  the  unsteady  aerodynamic  forces  acting  on  an 
aircraft.  As  a  result,  flutter  testing  of  a  scaled  model  in  a  transonic  wind 
tunnel  is  always  used  to  generate  corrections  to  flutter  speeds  computed  by  linear 
methods.  However,  the  design,  construction  and  testing  of  a  wind  tunnel  flutter 
model,  and  the  analysis  of  the  resulting  data,  require  over  a  year’s  time.  For 
this  reason,  leading  authorities  in  this  field  have  recently  suggested  [6]  that  “ The 
results  of  a  finite  number  of  [nonlinear]  CFD  [computational  fluid  dynamics] 
solutions  could  be  used  as  a  replacement  for  wind  tunnel  testing,  assuming  a 
validated  code  was  available.”  Coincidentally,  this  paper  addresses  this  very 
same  issue  in  the  context  of  an  F-16  Block  40  fighter. 

More  specifically,  our  objectives  here  are  to  overview  a  three-field  formulation 
of  fluid-structure  interaction  problems  for  nonlinear  computational  aeroelastic- 
ity,  assess  the  latest  advances  in  the  corresponding  CFD  and  CSM  (compu¬ 
tational  structural  mechanics)  algorithms,  briefly  describe  a  high-performance 
simulation  software  developed  at  the  University  of  Colorado  based  on  these 
ideas,  report  on  its  validation  with  flight  test  data  for  an  F-16  Block  40  con¬ 
figuration,  and  discuss  the  feasibility  as  well  as  the  merit  of  this  computational 
strategy  for  extracting  accurately  the  flutter  envelopes  of  high  performance  civil¬ 
ian  and  military  aircraft.  To  this  effect,  we  organize  the  remainder  of  this  paper 
as  follows. 

In  Section  2,  we  overview  a  three-field  formulation  of  coupled  fluid-structure 
interaction  problems  that  is  flexible  enough  to  accommodate  all  types  of  non- 
linearities.  In  Section  3,  we  identify  the  main  computational  'ssues  associated 
with  this  formulation.  For  some  of  these  issues,  we  take  a  neT"  look  at  the  lat¬ 
est  developments  and  achievements,  and  highlight  their  impact  on  the  overall 
performance  of  the  solution  of  nonlinear  aeroelastic  problems.  In  Section  4, 
we  briefly  describe  the  CFD/ CSM  based  simulation  software  developed  at  the 
University  of  Colorado  for  solving  general  aeroelastic  problems.  This  software 
is  based  on  the  three-field  formulation  and  incorporates  all  the  advances  in  the 
computational  methods  described  in  this  paper.  In  Section  5,  we  apply  this 
simulation  software  to  the  prediction  of  the  aeroelastic  frequencies  and  damp¬ 
ing  coefficients  of  an  F-16  Block  40  configuration  in  various  subsonic,  transonic, 
and  supersonic  regimes.  For  this  purpose,  we  represent  first  this  fighter  by  its 
typical  wing  section.  Then,  we  consider  a  three-dimensional  detailed  finite  ele¬ 
ment  (FE)  model  of  its  structure.  In  both  cases,  we  comment  on  the  obtained 
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numerical  results  and  compare  them  to  the  measured  flight  test  data.  Finally  in 
Section  6,  we  conclude  this  paper  by  a  discussion  of  the  feasibility  and  merit  of 
our  CFD/CSM  based  methodology  for  the  flutter  analysis  of  high  performance 
aircraft  in  the  transonic  and  other  regimes. 


2  Three-field  formulation  of  nonlinear  aeroelas- 
tic  problems 

First,  we  overview  a  three-field  formulation  of  coupled  fluid-structure  interac¬ 
tion  problems  that  was  first  introduced  in  [7].  This  formulation  is  quite  general. 
It  can  address  many  aeroelastic  problems  besides  flutter,  including  the  predic¬ 
tion  of  steady  and  unsteady  loads  as  well  as  control  surface  effects  in  level  flight 
as  well  as  during  maneuvering,  aeroelastic  tailoring,  and  performance  analysis. 
In  this  formulation,  the  structure  is  no  longer  restricted  to  a  harmonic  motion 
with  small  displacement  amplitudes,  and  is  not  necessarily  represented  by  a 
truncated  basis  of  its  normal  modes.  In  principle,  there  is  also  no  reason  to 
confine  its  constitutive  modeling  to  that  of  an  elastic  material.  However,  while 
aircraft  structures  can  undergo  large  displacements  and  rotations,  they  seldom 
experience  large  strains.  Therefore  in  many  applications,  the  nonlinear  modeling 
of  the  structural  behavior  can  be  limited  to  the  proper  accounting  of  nonlinear 
geometric  and  free  play  effects.  More  importantly,  the  aerodynamic  forces  act¬ 
ing  on  the  structure  are  no  longer  predicted  in  this  formulation  by  the  use  of  a 
linear  aerodynamic  operator,  because  of  the  limitations  associated  with  such  an 
approach,  particularly  in  the  transonic  regime.  Rather,  these  unsteady  forces 
axe  determined  from  the  solution  of  the  compressible  Euler  equations  when  vis¬ 
cous  effects  can  be  neglected,  and  the  solution  of  the  compressible  Navier-Stokes 
equations  augmented  by  a  large  eddy  simulation  or  turbulence  model  otherwise. 
Furthermore,  no  restriction  is  imposed  on  the  nature  of  the  fluid-structure  cou¬ 
pling,  which  is  numerically  modeled  by  suitable  fluid-structure  interface  bound¬ 
ary  (or  transmission)  conditions.  One  difficulty  in  handling  numerically  the 
fluid-structure  coupling  stems  from  the  fact  that  the  structural  equations  are 
usually  formulated  with  material  (Lagrangian)  coordinates,  while  the  fluid  equa¬ 
tions  are  typically  written  using  spatial  (Eulerian)  coordinates.  Therefore,  a 
straightforward  approach  to  the  solution  of  the  coupled  fluid-structure  dynamic 
equations  requires  moving  at  each  time-step  at  least  the  portions  of  the  fluid  grid 
that  axe  close  to  the  moving  and  flexing  aircraft.  This  can  be  appropriate  for 
small  displacements  of  the  structure  but  may  lead  to  severe  grid  distortions  when 
it  undergoes  large  motion.  Different  approaches  have  emerged  as  an  alternative 
to  partial  regridding  in  transient  aeroelastic  computations,  among  which  stand 
out  the  arbitrary  Lagrangian-Eulerian  (ALE)  formulation  [8]  and  the  closely 
related  method  of  dynamic  meshes  [9].  These  approaches  treat  a  computational 
aeroelasticity  problem  as  a  two-field  coupled  problem.  However,  a  moving  mesh 
(Fig.  1)  can  also  be  viewed  as  a  pseudo-structural  (or  fictitious  structural)  sys¬ 
tem  with  its  own  behavior  [7],  and  therefore,  the  coupled  transient  aeroelastic 
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problem  can  be  formulated  as  a  three-  rather  than  two-field  problem:  the  fluid, 
the  structure,  and  the  dynamic  fluid  mesh.  This  three-field  formulation  has  shed 
new  light  on  the  mathematical  understanding  of  the  numerical  behavior  of  var¬ 
ious  algorithms  applied  to  the  solution  of  the  coupled  fluid-structure  problem, 
and  has  enabled  the  development  of  faster  solution  algorithms  [10-14]. 

Clearly,  the  simultaneous  solution  of  the  governing  nonlinear  fluid,  fluid 
mesh,  and  structure  equations  of  motion  is  computationally  intensive,  and  raises 
some  concerns  about  the  feasibility  and  practiced  usefulness  of  this  approach  in 
production  environments.  An  important  objective  of  this  paper  is  to  show  that 
because  of  significant  advances  in  computational  methods  and  the  advent  of 
parallel  processing,  the  three-field  and  CFD/CSM  based  solution  of  aeroelas- 
tic  problems  is  now  sufficiently  mature  and  fast  to  be  considered  at  least  as 
a  reliable  simulation  environment,  if  not  to  be  adopted  as  a  design  tool,  for 
addressing  the  critical  flight  conditions  of  a  high  performance  aircraft. 

2.1  Governing  multidisciplinary  equations 

A  fluid-structure  interaction  problem  can  be  described  by  the  following  coupled 
partial  differential  equations 

JVx.R{w)  (3a) 

b  (3b) 

0  (3c) 

Equation  (3a)  is  the  ALE  conservative  form  of  the  Navier-Stokes  equations. 
Here,  t  denotes  the  time,  x(£)  denotes  the  time-dependent  position  or  displace¬ 
ment  of  a  fluid  grid  point  (depending  on  the  context  of  the  sentence  and  the 
equation),  f  its  position  in  a  reference  configuration,  J  =  det(dx/d £),  w  is  the 
fluid  state  vector  using  the  conservative  variables,  and  F  and  R  denote  respec¬ 
tively  the  convective  and  diffusive  ALE  fluxes.  Equation  (3b)  is  the  elastody- 
namic  equation  where  us  denotes  the  displacement  field  of  the  structure  and  ps 
its  density,  as  and  es  denote  respectively  the  stress  and  strain  tensors,  and  b 
represents  the  body  forces  acting  on  the  given  structure.  Equation  (3c)  governs 
the  dynamics  of  the  moving  fluid  grid.  It  is  similar  to  the  elastodynamic  equa¬ 
tion  because  the  dynamic  mesh  is  viewed  here  as  a  pseudo-structural  system.  A 
tilde  notation  is  used  to  designate  the  fictitious  mechanical  quantities  [1].  For 
the  sake  of  notational  simplicity,  the  various  Dirichlet  and  Neumann  boundary 
conditions  intrinsic  to  each  of  the  fluid  and  structure  problems  are  omitted. 

Equation  (3a)  and  equation  (3c)  are  directly  coupled.  If  up  denotes  the  ALE 
displacement  field  of  the  fluid  and  p  its  pressure  field,  <75  and  a p  the  structure 
stress  tensor  and  the  fluid  viscous  stress  tensor,  T  the  fluid-structure  interface 
boundary  (wet  boundary  of  the  structure),  and  n  the  normal  at  a  point  to  T, 
the  fluid  and  structure  equations  are  coupled  by  the  interface  boundary  —  or 


d(Jw ),  .  .  ox 

Ps-jfif  ~  dM<Ts(e5(us),  ^f(“s)))  = 


dt 

-  div(cr(l(x)))  = 
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transmission  —  conditions 


on  r 


(4a) 

(4b) 


as-Ti  =  —pn  +  <7F-n 


dus  _  duF 
dt  dt 


on  r 


The  first  of  these  two  transmission  conditions  states  that  the  tractions  on  the 
wet  surface  of  the  structure  axe  in  equilibrium  with  those  on  the  fluid  side  of  T . 
The  second  of  Eqs.  (4)  expresses  the  compatibility  between  the  velocity  fields 
of  the  structure  and  the  fluid  at  the  fluid-structure  interface.  For  inviscid  flows, 
this  second  equation  is  replaced  by  the  slip  wall  boundary  condition 


dup 

~df'n 


dus 

~dt'U 


on  r 


(5) 


The  equations  governing  the  structure  and  dynamic  mesh  motions  are  coupled 
by  the  continuity  conditions 


X  —  Us 

dx  _  dus 
~dt  ~  ~dt 


on  r 
on  r 


(6a) 

(6b) 


2.2  Semi-discretization  of  the  governing  equations 

The  spatial  approximation  of  the  ALE  conservative  form  of  the  Navier-Stokes 
equations  by  FE  or  finite  volume  (FV)  schemes  leads  to  semi-discrete  equations 
that  can  be  written  as 

(V(x)w)  +  F{w,x,x)  =  R(w,x)  (7) 

where  a  bold  font  designates  the  discrete  counterpart  of  a  field  variable,  and  a 
dot  a  time  derivative.  In  the  case  of  a  FV  semi-discretization,  V  is  the  matrix 
of  the  cell  volumes  and  F  and  R  are  respectively  the  numerical  convective  and 
diffusive  fluxes  approximating  the  integral  of  the  physical  flux  function  over  the 
cell  interfaces. 

The  semi-discretization  by  FE  of  the  structural  equations  of  dynamic  equi¬ 
librium  leads  to1 

Mils  +  fT(us,us)  =  fs(us,w)  +  ft  (8) 

where  M  is  the  FE  lumped  mass  matrix,  us  the  generalized  displacement  vec¬ 
tor,  f™1  the  vector  of  internal  forces,  fa$  the  vector  of  aerodynamic  forces 
applied  on  the  wet  surface  of  the  structure,  and  the  vector  of  other  exter¬ 
nal  forces  acting  on  the  structure.  If  for  the  problem  of  interest  the  structure 
remains  linear,  Eq.  (8)  simplifies  to 

Mus  +  Cus  +  Kus  =  fse  (us,  w)  +  / esxt  (9) 

1This  specific  expression  assumes  that  the  rotational  inertia  forces  are  negligible. 
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where  C  and  K  denote  the  FE  damping  and  stiffness  matrices. 

Let  the  subscript  i  designate  the  grid  points  located  in  the  interior  of  a 
computational  domain,  and  the  subscript  b  designate  those  grid  points  located 
at  the  fluid-structure  interface  T.  If  the  dynamic  mesh  is  assimilated  with  a 
quasi-static  pseudo-structural  system,  the  semi-discrete  equation  governing  the 
evolution  of  the  dynamic  mesh  can  be  written  as 

KuXi  =  -Kibxb  xb  =  Uush  (10) 

where  K  is  the  fictitious  time-depeudent  stiffness  matrix  resulting  from  the  FE 
semi-discretization  of  Eq.  (3c)  [1],  and  U  is  a  transfer  matrix.  If  the  fluid  and 
structure  meshes  have  compatible  interfaces,  U  =  I.  Otherwise,  U  is  given  by 
the  FE  discretization  of  the  second  of  the  interface  conditions  (4). 

3  Computational  issues  and  advances 

In  the  linear  theory  of  aeroelasticity,  the  air  surrounding  a  flying  aircraft  can  be 
interpreted  as  an  “algebraic”  damper  whose  sign  depends  on  the  flight  condi¬ 
tions.  When  positive,  it  attenuates  any  aircraft  vibration  excited  by  some  initial 
disturbance.  When  zero,  it  only  entertains  it,  and  when  negative,  it  amplifies 
that  vibration.  In  other  words,  depending  on  the  flight •  conditions  and  par¬ 
ticularly  the  Mach  number,  the  air  surrounding  a  vibrating  aircraft  can  either 
extract  energy  from  it,  or  act  as  a  neutral  agent  towards  it,  or  feed  it  energy  and 
cause  it  to  flutter.  This  energy  interpretation  of  the  flutter  mechanism  under¬ 
scores  the  importance  of  conserving  as  much  as  possible  the  energy  transferred 
between  the  fluid  and  structure  subsystems  when  discretizing  the  transmission 
conditions  (4)  and  solving  the  coupled  system  of  equations  (7,  8,  10).  Indeed, 
the  extraction  (transmission)  from  (to)  the  structure  across  the  fluid-structure 
interface  T  of  any  significant  amount  of  spurious  numerical  energy  can  artifi¬ 
cially  stabilize  (destabilize)  an  otherwise  unstable  (stable)  aeroelastic  system. 
Since  the  three-field  aeroelastic  problem  (7,  8,  10)  is  formulated  in  the  time  do¬ 
main,  extracting  the  wet  frequencies  and  damping  coefficients  of  the  underlying 
structure  for  flutter  analysis  require  post-processing  the  numerical  output  — 
for  example,  the  displacement  field  v  ■;(£)  —  by  a  parameter  identification  algo¬ 
rithm.  Computational  efficiency  suggests  using  for  that  purpose  an  identifica¬ 
tion  algorithm  that  requires  as  few  cycles  as  possible  of  the  predicted  structural 
response.  This  in  turn  underscores  the  importance  of  producing  a  sufficiently 
accurate  short  window  of  the  time-response,  and  therefore  time-integrating  the 
coupled  fluid-structure  and  not  individual  fluid  and  structure  equations  of  mo¬ 
tion  (7,  8)  with  a  second-order  time-accurate  scheme.  Usually,  the  aeroelastic 
response  of  the  structure  is  dominated  by  its  low  modes.  For  this  reason,  com¬ 
putational  speed  —  which  is  essential  for  production  environments  —  favors 
implicit  schemes  and  large  computational  steps,  which  underscores  the  impor¬ 
tance  of  paying  special  attention  to  the  numerical  stability  properties  of  the 
scheme  designed  for  time-integrating  the  coupled  fluid-structure  equations  of 
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motion.  Next,  we  review  some  recent  achievements  in  these  areas.  Because 
of  space  limitation,  we  focus  on  discussing  results,  and  refer  the  reader  to  the 
appropriate  literature  for  details. 

3.1  CFD  on  moving  unstructured  grids 

The  governing  fluid  equation  (7)  differs  from  the  standard  FE  or  FV  semi- 
discretization  of  the  Navier-Stokes  equations  in  that  it  is  formulated  on  a  mov¬ 
ing  rather  than  a  fixed  grid.  Therefore,  the  time-discretization  of  this  equa- 

tion,  which  incurs  the  approximation  of  the  integrals  ftn  £‘(u i,x,x)dt  and 

ftn  R(w,x)  dt,  raises  the  question  of  where  to  evaluate  the  convective  and 
diffusive  fluxes  when  the  grid  moves  from  its  position  xn  at  time  tn  to  a  posi¬ 
tion  xn+1  at  time  £n+1 .  A  straightforward  answer  to  this  question  is  to  evaluate 
these  fluxes  on  the  mesh  configuration  xn  when  the  chosen  time-integrator  is 
explicit,  and  xn+1  when  the  time-integrator  is  implicit.  For  small  time-steps 
Atn  =  £n+1  -  tn,  it  may  not  matter  where  the  fluxes  axe  evaluated  because 
in  that  case  the  difference  between  the  mesh  configurations  xn  and  xn+l  may 
be  insignificant.  However,  for  the  large  time-steps  dictated  by  computational 
efficiency,  the  method  of  evaluation  of  the  integrals  Jtn  F(w,x,x)  dt  and 

ffn  R(w,x)dt  can  have  a  dramatic  effect  on  the  performance  of  the  time- 
integration  of  the  governing  fluid  equation  (7).  To  address  this  issue,  it  was 
proposed  in  [10-13]  to  first  select  a  time-integrator  that  performs  well  on  fixed 
grids,  and  then  extend  it  to  moving  grids  by  evaluating  a  moving  flux  as  the 
time-average  of  a  certain  number  of  fluxes  computed  on  a  suite  of  carefully 
chosen  mesh  configurations.  For  example,  the  classical  three-point  backward 
difference  scheme  fits  well  the  selection  criteria  stated  in  the  introduction  of 
this  section:  it  is  implicit,  second-order  time-accurate  on  fixed  grids,  and  un¬ 
conditionally  stable  for  the  usual  test  problem  on  fixed  grids.  Its  extension  to 
moving  grids  as  proposed  in  [10-13]  can  therefore  be  written  as  follows 

an+iVn+1wn+1  +anVnwn  +  an-1Vn~1wn~1 

n°, 

+  AtnYl^F(wn+\x‘,xl)-AtnY/ciR(w'-+1,xds)=0  (11) 

5=1  5=1 


where 


C*n+1 


1  +  2  r 

1  +  T  ’ 


an  —  1  T,  Qn-1  — 


1  +  r’ 


A  tn 
At"-1 


and  V"  =  V (xn).  ncs  and  nd  axe  two  small  integers,  while  c£  and  cd  are  real 
coefficients  (averaging  weights)  that  satisfy  £  ccs  =  1  and  £  cd  =  1.  xcs,  xds, 

5=1  5=1 

and  xc3  are  linear  combinations  of  the  mesh  configurations  { xn~l , ..,  xn, ..,  xn+m  } 
and  velocities  {in“J, ±n, xn+q }  for  some  given  integers  Z,  m,  j,  and  q.  Note 
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that  Eq.  (11)  anticipates  two  different  time-averaging  procedures  for  the  moving 
convective  and  viscous  fluxes. 

The  computational  complexity  of  the  scheme  (11)  can  be  reduced  by  aver¬ 
aging  the  mesh  configurations  themselves  instead  of  the  fluxes  associated  with 
them,  which  leads  to 

a„+i  Vn+1mn+1  +  anVnwn  +  an-1Vn-1ton-1 

n‘  nc,  nf 

+  AtnF(wn+\^clxl,^clxl)-AtnR(wn+\Y.cdsxds)  =  0  (12) 

5=1  5=1  S=1 

This  alternative  approach  was  proposed  and  discussed  in  [13]  in  the  context 
of  the  FV  semi-discretization  of  the  governing  flow  equations.  It  requires  the 
computation  of  a  single  convective  flux  and  a  single  diffusive  flux  per  time-step, 
whereas  the  approach  summarized  in  Eq.  (11)  requires  the  computation  at  each 
time- step  of  ncs  convective  and  nd  diffusive  fluxes.  In  the  remainder  of  this 
paper,  we  adopt  the  form  (12)  of  the  extension  to  moving  grids  of  the  three- 
point  backward  difference  scheme.  However,  most  of  the  results  we  present  for 
that  form  also  hold  for  the  version  (11)  of  that  extension. 

Whether  scheme  (11)  or  scheme  (12)  is  chosen  for  extending  the  three-point 
backward  difference  scheme  to  moving  grids,  it  remains  to  determine  ncs  and  nf , 
the  averaging  coefficients  ccs  and  cd ,  the  mesh  configurations  xcs ,  and  xd,  and 
the  mesh  velocities  xc3 .  To  this  effect,  two  approaches  can  be  adopted.  The  first 
one  consists  in  choosing  these  parameters  so  that  the  resulting  time-integrator 
satisfies  its  corresponding  discrete  geometric  conservation  law  (DGCL)  [15,16]. 
The  second  approach  consists  in  selecting  these  parameters  so  that  the  resulting 
time-integrator  remains  second-order  time-accurate  on  moving  grids. 

The  so-called  DGCL  is  a  law  which  states  that  the  computation  of  the  geo¬ 
metric  parameters  associated  with  a  moving  grid  should  be  performed  in  such  a 
way  that,  independently  of  the  mesh  motion,  the  numerical  scheme  constructed 
for  time-integrating  the  flow  equations  on  a  moving  grid  preserves  the  state  of 
a  uniform  flow.  The  idea  of  computing  the  discrete  mesh  velocities  and  other 
geometric  parameters  as  to  preserve  a  certain  physical  quantity  goes  back  to  the 
early  days  of  CFD.  The  terminology  “geometric  conservation  law”  was  coined  in 
1979  by  Thomas  and  Lombard  [17]  who  derived  this  concept  from  the  law  of  mass 
conservation  in  a  spatial  region  bounded  by  a  moving  surface,  and  applied  it  to 
construct  an  improved  finite  difference  method  for  flow  computations  on  mov¬ 
ing  grids.  This  concept  was  subsequently  extended  to  characterize  geometrically 
conservative  schemes  as  algorithms  that  preserve  the  entire  state  of  a  uniform 
flow.  First-order  time-accurate  and  geometrically  conservative  ALE  FV  schemes 
were  presented  and  discussed  in  [11,12,18].  First-order  time-accurate  and  geo¬ 
metrically  conservative  ALE  FE  schemes  were  presented  in  [11,12].  DGCLs  for 
second-order  time-accurate  ALE  FV  schemes  have  also  been  developed  and  dis¬ 
cussed  in  [13],  For  example,  parameterizing  the  sought-after  mesh  configuration 
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(13a) 

(13b) 


xcs  and  its  corresponding  velocity  field  as  suggested  in  [13] 

**  =  tf+1Xn+1  +  tfx"  +  (1  -  t£+1  -  J7?)xfl-1 

±;  =  e?+1xn+1  +  e*xn  -  (0sn+1  +  ens)xn~l 


and  requiring  that  the  scheme  (12)  satisfies  its  DGCL  leads  to  a  nonlinear 
system  of  equations  in  the  unknowns  nc„  ccs,  {j?"+p}£=o>  {^"+p}p=o-  bi  ^e 

case  of  three  dimensions  and  a  FV  semi-discretization  method,  one  of  the  many 
possible  solutions  of  that  nonlinear  system  of  equations  leads  to  (see  [13]) 


nc. 


=  4;  c\ 

1(1 


■<% 


Qn^l  . 
2 


{x 


Cg  —  ^  — 

^  =  |(l-^n+1  +  |(l  +  ^n 
*1  =  Id  +  ^)*n+1  +  Id  -  Ts)*" 
^ = 1(1  -  + 1(1 + 

*c4  =  Id  +  Tf)*"  +  5(1  -  Ts)*”-1 


(14) 


This  completes  the  description  of  one  instance  of  the  scheme  (12)  that  is  based  on 
a  four-point  discretization  of  the  moving  convective  fluxes.  Enforcing  a  DGCL 
cannot  determine  the  coefficients  cd  and  mesh  configurations  xd,  because  the 
moving  viscous  fluxes  vanish  for  a  uniform  flow.  To  determine  these  parameters, 
one  can  adopt  the  second  approach  mentioned  above  and  described  next. 

Alternatively,  the  parameterization  (13)  and  a  similar  parameterization  of 
the  sought-after  mesh  configuration  x d  can  be  substituted  into  the  scheme  (12). 
Next,  the  local  truncation  error  of  the  resulting  time-integration  algorithm  can 
be  computed  as  usual  by  applying  Taylor’s  expansion  around  the  variables  at 
tn.  Finally,  the  unknowns  c£,  cdy  and  the  parameters  governing  the  mesh  con¬ 
figurations  xc$ ,  xcs  and  xd  can  be  determined  by  requiring  that  the  zero-,  first-, 
and  second-order  terms  of  the  local  truncation  error  vanish.  It  turns  out  that 
this  second  approach  for  determining  c8,  cd,  x xcsy  and  xds  works  equally  well 
when  ncs  is  set  a  priori  to  nc8  =  4,  nd  and  cd  are  set  a  priori  to  nd  =  1  and  cd  =  1, 
and  xd  is  parameterized  differently  than  xca  using  the  simpler  linear  map 

xd3  =  e+1zn+1  +  (1  —  C+1)*n  (15) 

which  was  first  proposed  in  [22].  In  that  case,  requiring  that  the  scheme  (12) 
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be  second-order  time-accurate  on  moving  grids  leads  to 


4 


(16a) 

S=1 

5=1 

(16b) 

JL  1  —  rjn+1  —  7)n 

i>sw+i  -  \  jk) = i 

S=1 

(16c) 

5— 1 

(16d) 

C+1  =  1 

(16e) 

Again,  the  system  of  equations  (16)  admits  multiple  solutions,  among  which 
we  note  two  solutions  corresponding  to  the  following  averaging  coefficients  and 
mesh  configurations 


rn%  =  4;  c\  =  1;  c§  =  c§  =  c%  =  0 
x\  =  xn+1 

,  — c  _  an+1Xn+1+aTiXn+an-ia:n-1 
\  Xl  -  At" 


(17) 


and 


Ctn-l 

2r 


(18) 


It  follows  that  at  least  two  different  instances  of  the  parameterized  time-integrator 
(12)  achieve  second-order  accuracy  on  moving  grids.  The  first  one  is  based  on 
the  one-point  rule  described  by  Eqs.  (17).  This  algorithm  does  not  however 
satisfy  its  DGCL  (for  a  proof,  the  reader  can  simply  verify  that  Eqs.  (17)  do 
not  satisfy  the  DGCL  requirements  of  the  time-integrator  (12)  which  are  listed 
in  [13]).  The  second  instance  of  the  parameterized  time-integrator  (12)  is  based 
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on  the  four-point  rule  described  by  Eqs.  (18),  which  is  identical  to  that  de¬ 
scribed  by  Eqs.  (14);  hence,  this  second  instance  of  scheme  (12)  satisfies  its 
DGCL.  Therefore,  using  what  we  have  referred  to  as  the  second  approach  for 
extending  the  three-point  backward  difference  scheme  to  moving  grids,  one  can 
construct  different  time-integrators  for  the  governing  semi-discrete  fluid  equa¬ 
tions  (7)  that  are  all  second-order  time-accurate,  but  which  do  not  necessarily 
satisfy  their  DGCL.  This  is  consistent  with  the  following  theoretical  result  es¬ 
tablished  by  Guillard  and  Farhat  [15,16]:  “for  a  given  scheme  that  is  p-order 
time-accurate  on  a  fixed  mesh,  satisfying  the  corresponding  p-discrete  geometric 
conservation  law  is  [only]  a  sufficient  condition  for  this  scheme  to  be  at  least 
first-order  time-accurate  on  a  moving  mesh”.  Therefore,  the  question  becomes 
whether  there  is  any  particular  reason  to  discriminate  between  those  second- 
order  extensions  to  moving  grids  of  the  three-point  backward  difference  scheme 
that  satisfy  their  DGCLs,  and  those  that  do  not  satisfy  it? 

In  [10,12],  it  was  shown  numerically  that  for  typical  aeroelastic  computa¬ 
tions,  violating  the  DGCL  can  introduce  a  parasitic  weak  instability  in  the  lift 
response,  and  can  lead  to  the  prediction  of  an  erroneous  flutter  speed.  Motivated 
by  these  observations,  Formaggia  and  Nobile  have  investigated  the  solution  of 
linear  advection-diffusion  problems  on  moving  grids  by  ALE  FE  methods  [20]. 
They  have  shown  that  for  these  linear  problems,  satisfying  the  corresponding 
first-order  discrete  geometric  conservation  law  is  a  sufficient  condition  for  the 
backward  Euler  implicit  scheme  to  be  unconditionally  stable.  This  new  result 
sheds  some  light  on  the  relationship  between  the  DGCL  and  numerical  stability. 
However,  it  does  not  take  into  account  the  nonlinearities  that  characterize  Eu¬ 
ler  flows,  and  stops  short  from  predicting  the  behavior  of  an  ALE  scheme  that 
does  not  satisfy  its  corresponding  DGCL.  For  this  reason,  Farhat,  Geuzaine 
and  Grandmont  have  investigated  further  the  theoretical  status  of  the  DGCL, 
and  exposed  its  relation  to  nonlinear  stability.  More  specifically,  using  a  three- 
dimensional  nonlinear  scalar  hyperbolic  conservation  law  as  a  model  problem, 
they  have  proved  for  sample  arbitrary  Lagrangian  Eulerian  schemes  that  the 
DGCL  requirement  corresponds  to  a  necessary  and  sufficient  condition  for  a 
numerical  scheme  to  achieve  nonlinear  stability.  Consequently,  they  have  also 
shown  that  an  ALE  scheme  which  violates  its  DGCL  is  bound  to  exhibit  spuri¬ 
ous  oscillations  and  overshoots  fcr  practical  computational  time-steps,  and  can 
occasionally  exhibit  an  unbounded  behavior.  For  all  these  reasons,  and  because 
the  computational  overhead  associated  with  enforcing  a  DGCL  is  minimal  —  for 
example,  whether  equipped  with  Eqs.  (17)  or  Eqs.  (18),  scheme  (12)  computes 
only  one  convective  and  one  diffusive  flux  per  time-step  —  we  prefer  numerical 
methods  that  satisfy  their  DGCLs.  In  particular,  we  prefer  the  four-point  rule 
version  (18)  of  the  time-integrator  (12)  over  its  one-point  rule  version  (17).  Fur¬ 
thermore,  we  point  out  that  in  [16]  and  references  cited  therein,  it  was  shown 
that  for  a  specified  accuracy,  the  four-point  rule  (18)  allows  a  computational 
time-step  that  is  more  than  an  order  of  magnitude  larger  than  that  afforded  by 
straightforward  and  geometrically  non  conservative  extensions  to  moving  grids 
of  the  three-point  backward  difference  scheme.  Hence,  it  is  now  well  established 
that  satisfying  the  governing  DGCL  improves  significantly  both  robustness  and 
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performance. 


3.2  Energy  conservative  exchange  of  aerodynamic  and  elas- 
todynamic  data 

Next,  we  turn  our  attention  to  the  discretization  of  the  transmission  conditions 
(4,6).  In  order  to  address  the  practiced  case  where  the  fluid  and  structure 
computational  domains  have  non-matching  discrete  interfaces,  we  denote  by  IV 
and  Ts  the  discrete  representations  of  the  fluid-structure  interface  T  on  the  fluid 
and  structure  sides,  respectively. 

The  energy  transferred  from  the  fluid  to  the  structure  through  T  during  the 
time-interval  [tn,  tn+1]  is 

rtn+1 

SEp71*1  =  —  J  ( J  (-pn  +  <7jr.n)x  ds'j  dt  (19) 

and  that  received  by  the  structure  during  the  same  time-interval  is 

t«+i 

Jf5n+1  =  J  (y  as. nits  da)  dt  (20) 

From  the  above  expressions  and  Eqs.  (4,6),  it  follows  that  the  discretization  of 
the  transmission  conditions  (4,6)  has  a  direct  impact  on  the  conservation  of  the 
energy  transferred  between  the  fluid  and  the  structure  through  T. 

Whichever  approximation  method  (interpolation,  projection,  etc.)  is  chosen 
for  enforcing  on  T  the  compatibility  of  the  displacement  fields  of  the  fluid  mesh 
and  the  structure  (6),  its  outcome  can  be  written  as  follows 

*s 

Xj  =  Cji  Ug.  j  €  IV,  i  €  Ts  (21) 

i=  1 

where  Xj  is  the  discrete  value  of  x  at  the  fluid  point  j,  and  us{  is  the  discrete 
value  of  us  at  the  structure  node  i .  The  integer  is  and  real  coefficients  Cji 
depend  on  the  chosen  method  of  approximation.  The  superscript  P  is  used  to 
designate  a  “prediction”  of  the  motion  of  the  structure.  If  the  huid  and  structure 
subsystems  are  solved  simultaneously,  tif  =  us.  On  the  o hand,  if  they  are 
advanced  in  time  by  a  staggered  procedure  (see  Section  3.3),  will  in  general 
differ  form  us{  because  of  a  time-lag  between  the  fluid  and  structure  partitions. 

Consider  now  a  virtual  displacement  field  x  that  is  zero  on  each  degree 
of  freedom  of  the  moving  fluid  grid  except  on  those  lying  on  the  boundary 
I>.  Whichever  method  is  chosen  for  approximating  Eq.  (3c)  and  therefore 
constructing  the  pseudo-structural  stiffness  matrix  K  (see  Eq.  (10)),  x  can  be 
expressed  as  follows 


3F 

x  —  ^  DjXj  j^TF 
3=1 


(22) 
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where  Dj  is  some  function  with  a  local  or  global  support  on  Tf.  From  Eq. 
(19)  and  Eq.  (22),  it  follows  that  the  virtual  work  during  [tn,  tn+1]  of  the  fluid 
tractions  acting  on  I>  is 


^n-f  1  j  jr 

6W£+1  =  £  (^  /  +  aF-n)Djij  ds)  dt 

.n  +  i  jF 
=  dt 
Jt"  ,=i 

where  $ j  has  the  physical  meaning  of  a  numerical  pressure  flux  and  is  given  by 

=  [  (-pn  +  aF-n)Dj  ds  (24) 

Jr* 


(23) 


Substituting  Eq.  (21)  into  Eq.  (23)  gives 

jn  +  l  jF 


is  rtn+1 


<5^f+1  =  [  £  Mi  dt  =  J2  [  /f.  «s, 

•'*"  i=i  <=i Jt" 


(25) 


where 


3F 


f  Fi  =  j  e  IV,  i  e  Ts 

j=l 


(26) 


If  the  nodal  aerodynamic  forces  (and  moments)  acting  on  a  structure  node  i 
lying  on  Ts  axe  denoted  by  the  virtual  work  during  [ tn ,  £n+1]  of  these 
forces  is 


is  is 

8WS+1  =  I  £  ft  us,  dt  =  J2  ts,  us,  dt 

Jtn  i=  1  i=l  Jtn 


(27) 


Conserving  the  transfer  of  energy  through  T  requires  enforcing  5W£+l  =  SWg*1 
for  any  pair  of  virtual  displacement  vectors  x  and  us  satisfying  (21).  From  Eqs. 
(25,27),  this  implies  enforcing 


/ 

Jtn 


f  +  1 

fsAsi  dt  =  I  fFAs> 
Jtn 


dt 


(28) 


The  evaluation  of  the  time-integrals  in  Eq.  (28)  depends  on  the  time-integration 
schemes  chosen  for  the  structure  and  fluid  analyses,  and  is  discussed  in  great 
details  in  [14].  Here,  we  focus  on  the  case  where  Eq.  (8)  is  time-integrated 
with  the  popular  midpoint  rule  —  which  is  second-order  time-accurate  and  in 
the  linear  case  identical  to  the  Newmark  algorithm  with  0  =  1/4  and  7  =  1/2 
[21]  —  and  the  flow  equations  (7)  are  time-integrated  with  the  second-order 
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time-accurate  extension  to  moving  grids  of  the  three-point  backward  difference 
scheme  described  by  Eqs.  (12,18).  In  that  case,  Eq.  (28)  becomes 


+  /!f 


(u 


n+1 

Si 


-«<?.) 


P 


(29) 


p 

If  the  fluid  and  structure  subsystems  axe  solved  simultaneously, 
and  Eq.  (29)  simplifies  to 


ft  =  fFi  =  f>;]M  J  erF,i&Ts  (30) 

J= 1 

In  that  case,  Eq.  (30)  describes  a  conservative  algorithm  for  transmitting  the 
fluid  forces  to  the  structure.  This  algorithm  is  independent  of  the  method  chosen 
for  discretizing  the  structure.  The  term  in  the  first  bracket  in  Eq.  (30)  depends 
exclusively  on  the  method  chosen  for  discretizing  the  flow  problem,  and  the  term 
in  the  second  bracket  depends  only  on  the  method  selected  for  transmitting  the 
displacement  of  the  structure  to  the  fluid  mesh. 

If  a  time-lag  is  introduced  in  the  solution  of  the  fluid  and  structure  sub¬ 
systems,  ug  tig.  In  that  case,  conserving  the  energy  transfer  through  T 

p  p 

requires  subiterating  on  ug  until  ug  =  ug.  However,  even  if  no  subiterations 
are  performed,  we  still  recommend  exchanging  the  aerodynamic  data  between 
the  fluid  and  the  structure  according  to  Eq.  (30),  and  discuss  in  Section  3.3  ef¬ 
fective  means  for  compensating  the  strict  loss  of  conservation  of  energy  transfer 
through  T. 


3.3  Higher-order  loosely  coupled  fluid-structure  time-integrators 

For  any  reasonably  detailed  FE  representation  of  the  aircraft  structure,  the 
simultaneous  solution  of  Eqs.  (3)  by  a  monolithic  scheme  is  neither  practical 
nor  software- wise  manageable.  For  this  and  other  reasons  related  to  compu¬ 
tational  efficiency,  a  partitioned  procedure  is  typically  employed  for  solving 
coupled  field  nonlinear  aeroelastic  problems.  In  such  a  procedure,  the  fluid 
and  structure  are  time-integrated  by  different  schemes  tailored  to  their  differ¬ 
ent  mathematical  models,  and  the  resulting  discrete  equations  are  solved  by 
a  staggered  algorithm  [10,23,24].  Such  a  strategy  simplifies  explicit /implicit 
treatment,  subcycling,  load  balancing,  software  modularity,  and  replacements 
as  better  mathematical  models  and  methods  emerge  in  the  fluid  and  structure 
disciplines.  The  most  popular  partitioned  procedure,  which  is  referred  to  in 
this  paper  as  the  Conventional  Serial  Staggered  (CSS)  procedure,  is  graphically 
depicted  in  Fig.  2.  Its  generic  cycle  between  tn  and  £n+1  can  be  described  as 
follows:  (1)  predict  ug^1  =  ug,  transfer  the  corresponding  motion  of  the  wet 
boundary  of  the  structure  to  the  fluid  subsystem,  and  update  the  position  of  the 
moving  fluid  mesh  accordingly,  (2)  advance  the  fluid  subsystem  to  tn+1  using  a 
given  flow  time-integrator  and  compute  a  new  set  of  aerodynamic  forces  /£+1 
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acting  on  I>,  (3)  transfer  /sen+1  =  /£+1  to  the  structure  subsystem,  (4)  apply 
fagnJr  1  to  the  structure  and  advance  it  to*  £n+1  using  a  given  time-integrator. 
Such  a  staggered  procedure,  which  can  be  described  as  a  loosely  coupled  solu¬ 
tion  algorithm,  can  also  be  equipped  with  a  subcycling  strategy  where  the  fluid 
and  structure  subsystems  are  advanced  using  different  time-steps  A tp  and  A t$- 

Unfortunately,  the  time-accuracy  of  the  CSS  procedure  is  in  general  at  least 
one  order  lower  than  that  of  its  underlying  flow  and  structure  time- integrators, 
and  its  numerical  stability  is  more  restrictive  than  that  of  the  flow  and  structure 
solvers.  Consequently,  its  maximum  allowable  time-step  is  much  smaller  than 
required  for  accuracy  purposes,  which  makes  it  a  slow  algorithm.  To  improve  the 
performance  of  this  simple  partitioned  procedure,  several  ad-hoc  strategies  have 
been  proposed  in  the  literature.  Essentially,  these  strategies  insert  some  type  of 
predictor/corrector  iterations  within  each  cycle  of  the  computations  in  order  to 
compensate  for  the  time-lag  between  the  fluid  and  structure  solvers  [23,25,26]. 
These  iterations  help  closing  the  gap  between  u£+1  and  ti5+1,  and  therefore 
as  explained  in  Section  3.2  (see  Eq.  (29))  they  also  help  conserving  the  transfer 
of  energy  between  the  fluid  and  structure  through  T.  As  a  result,  they  increase 
the  maximum  allowable  time-step  of  the  CSS  procedure.  However,  since  the 
numerical  complexity  of  each  predictor/corrector  iteration  is  almost  the  same 
as  that  of  one  cycle  of  staggering,  little  CPU  is  saved  by  such  enhancement 
strategies. 

Recently,  Piperno  and  Farhat  [14]  have  presented  an  alternative  approach 
for  improving  the  maximum  allowable  time-step  of  the  CSS  procedure  that  does 
not  increase  its  computational  cost  per  cycle.  Their  approach,  which  assumes 
that  /£+1  is  computed  as  in  Eq.  (26),  is  based  on  introducing  two  compu¬ 
tationally  economical  factors  for  compensating  the  time-lag  between  the  fluid 
and  the  structure  subsystems:  (1)  a  non-trivial  prediction  ^  and 

(2)  a  non  necessarily  trivial  transfer  of  the  aerodynamic  forces  to  the  struc¬ 
ture  fg6  /f+1*  More  specifically,  they  have  shown  that  given  two  time- 
integration  schemes  for  the  fluid  and  structure  equations  of  motion,  u^+l  and 
/r+1  can  be  designed  to  achieve 

NT/ At 

5W£  +  SWS  =  O  (A tp)  (31) 

n=0 

where  N  is  an  arbitrary  integer,  T  is  the  period  of  an  assumed  harmonic  vibra¬ 
tion  of  the  structure,  and  At  =  At?  =  A ts  is  a  fixed  time-step  for  the  fluid  and 
structure  time-integrators.  In  other  words,  ti£+1  and  fa$  can  be  designed 
to  construct  a  p-order  “energy-transfer-accurate”  CSS  procedure.  The  higher  p 
is,  the  closer  is  the  CSS  procedure  to  conserving  the  transfer  of  energy  through 
the  fluid-structure  interface.  For  example,  consider  the  case  where  the  CSS 
procedure  is  equipped  with  the  midpoint  rule  for  time-integrating  the  structure 
subsystem,  and  the  extension  to  moving  grids  of  the  three-point  backward  dif¬ 
ference  scheme  described  by  Eqs.  (12,18)  for  integrating  the  fluid  subsystem. 
When  u^lP  =  wg  and  fa£  +  =  /  j?+1i  this  CSS  procedure  is  only  first-order 
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energy- transfer- accurate.  On  the  other  hand,  if  is  predicted  by  the  fol¬ 

lowing  first-order  scheme 

uS+1'=»5  +  ^*S  (32) 

and  the  following  “improved”  vector  of  aerodynamic  forces  is  applied  on  the  wet 
surface  of  the  structure 

r/+1 = 2/r1  -  ff  (33) 

the  CSS  procedure  becomes  third  order  energy-transfer-accurate  and  its  max¬ 
imum  allowable  time-step  is  increased  by  a  factor  equal  to  5.  This  significant 
improvement  in  performance  is  achieved  without  any  predictor/corrector  itera¬ 
tion  and  almost  at  zero  cost. 

In  [14],  it  is  also  shown  that  the  Improved  Serial  Staggered  (ISS)  procedure 
proposed  by  Lesoinne  and  Farhat  in  [27]  is  also  third-order  energy-transfer- 
accurate,  and  allows  a  maximum  time-step  that  is  at  least  as  large  as  that  of  the 
CSS  procedure  equipped  with  the  first-order  prediction  (32)  and  the  improved 
aerodynamic  forces  (33).  This  is  not  surprising  because  as  shown  in  Fig.  3, 
the  ISS  procedure  incorporates  the  effects  of  the  improved  prediction  (32)  and 
aerodynamic  forces  (33). 

Finally,  we  note  that  the  combination  of  the  advanced  fluid  time-integrators 
highlighted  in  Section  3.1  and  improved  staggered  procedures  overviewed  herein 
has  improved  the  computational  efficiency  and  performance  of  the  state-of-the- 
art  technology  for  the  solution  of  nonlinear  transient  aeroelastic  problems  by  a 
factor  ranging  between  10  and  22  [27]. 


4  The  AERO-F/AERO-S  simulation  platform 

The  AERO-F,  AERO-S,  and  MATCHER  codes  developed  at  the  University  of 
Colorado  are  a  suite  of  software  modules  based  on  the  three-field  formulation 
described  in  this  paper  for  the  solution  on  nonlinear  transient  aeroelastic  prob¬ 
lems.  They  are  portable,  and  run  on  a  large  variety  of  computing  platforms 
ranging  from  Unix  workstations  to  shared  as  well  as  distributed  memory  mas¬ 
sively  parallel  computers. 

The  two-  and  three-dimensional  AERO-F  modules  model  a  flow  either  by  the 
Euler  equations,  or  the  averaged  Navier-Stokes  equations  equipped  with  the  k-e 
turbulence  model  and  a  wall  function.  They  operate  on  static  and  dynamic  un¬ 
structured  meshes.  More  specifically,  they  combine  a  Galerkin  centered  approxi¬ 
mation  for  the  viscous  terms,  and  a  Roe  upwind  scheme  for  the  convective  fluxes. 
Higher-order  spatial  accuracy  is  achieved  through  the  use  of  a  multidimensional 
piecewise  linear  reconstruction  that  follows  the  principle  of  the  Monotonic  Up¬ 
wind  Scheme  for  Conservative  Laws  [28].  Time-integration  on  fixed  grids  can 
be  performed  either  by  a  3-step  variant  of  the  explicit  Runge-Kutta  algorithm, 
or  by  the  implicit  three-point  backward  difference  scheme.  Time-integration  on 
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moving  grids  is  carried  out  as  described  in  this  paper.  All  linearized  systems  of 
equations  are  solved  by  the  RAS  (Restricted  Additive  Schwarz)  preconditioned 
GMRES  iterative  algorithm  [31]. 

The  AERO-F  modules  support  two  robust  structure  analogy  methods  for 
constructing  dynamic  meshes.  The  first  one  is  based  on  time-dependent  tor¬ 
sional  springs  [29,30].  The  second  method  is  based  on  the  total  Lagrangian 
approach  for  solving  a  fictitious  nonlinear  elasticity  problem  [32].  Both  meth¬ 
ods  share  in  common  the  idea  of  constructing  the  fictitious  stiffness  of  each  fluid 
mesh  element  so  that  it  increases  to  infinity  when  its  area  or  volume  decreases  to 
zero.  This  prevents  all  collapsing  mechanisms  (node-to-nodoy  node-to-edge,  and 
node-to-face)  from  occurring  during  the  mesh  motion.  For  applications  where 
the  structure  undergoes  large  rotations  —  for  example,  aircraft  maneuvering  — 
the  AERO-F  modules  invoke  a  corotational  scheme  to  accelerate  the  update  of 
the  mesh  motion  [33]. 

The  AERO-S  suite  of  structural  and  thermal  modules  are  capable  of  linear 
and  geometrically  nonlinear  static,  sensitivity,  vibration  (eigen),  and  transient 
FE  analyses  of  restrained  as  well  as  unrestrained  homogeneous  and  composite 
structures. 

The  AERO-F  and  AERO-S  suite  of  codes  communicate  via  run-time  soft¬ 
ware  channels.  They  exchange  aerodynamic  and  aeroelastic  data  across  non¬ 
matching  fluid  and  structure  mesh  interfaces  as  described  in  this  paper.  For  that 
purpose,  they  are  guided  by  information  generated  in  a  preprocessing  phase  by 
the  MATCHER  software  [34]. 

5  Application  to  the  prediction  of  the  aeroelas¬ 
tic  parameters  of  an  F-16  configuration 

The  validation  of  the  AERO-F,  AERO-S,  and  MATCHER  codes  for  the  flutter 
analysis  of  the  AGARD  Wing  445.6  is  described  in  [27].  For  this  relatively  simple 
wing  problem,  these  codes  and  the  underlying  algorithms  described  in  this  paper 
have  proved  to  be  capable  of  capturing  correctly  the  transonic  dip.  They  have 
also  demonstrated  a  superior  computational  efficiency  by  operating  accurately 
with  a  fluid  time-step  A tp  that  is  10  to  22  times  larger,  and  ^coupling  time-step 
At  =  max(AtS)  At^)  that  is  20  to  46  times  larger,  than  repo^d  in  the  literature 
for  other  nonlinear  aeroelastic  codes  [27,35].  Here,  we  validate  our  nonlinear 
transient  aeroelastic  simulation  technology  for  an  industrial  problem  in  view  of 
assessing  its  potential  for  replacing  flutter  testing  of  scaled  models  in  transonic 
wind  tunnels.  For  this  purpose,  we  simulate  the  flutter  clearance  of  an  F-16 
Block-40  in  clean  wing  configuration  but  with  tip  missiles,  for  0.8  <  Moo  <  1.4 
and  at  the  altitude  of  3,000  m.  For  this  purpose,  we  first  represent  the  structural 
behavior  of  this  fighter  by  that  of  its  typical  wing  section,  then  by  that  of  a 
detailed  three-dimensional  FE  model.  We  perform  all  computations  in  double 
precision  arithmetic  on  a  Silicon  Graphics  Origin  2000  parallel  computer.  We 
contrast  the  results  obtained  using  both  approaches  and  compare  them  to  flight 
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test  data  provided  by  the  Flight  Test  Center  at  the  Edwards  Air  Force  Base. 

5.1  Typical  wing  section  model 

From  a  structural  viewpoint,  a  typical  wing  section  (TWS)  is  characterized  by 
two  degrees  of  freedom  (dof)  h  and  9  which  represent  the  behavior  of  the  wing 
in  bending  and  torsion,  respectively  (see  Fig.  4).  It  is  identified  by  the  following 
parameters 


Ptws  =  {m,  h,%G,  {xg  -  %c),  Khy  Kq } 

where  m  is  the  total  mass  per  unit  span  of  the  wing,  xq  and  xc  designate  the 
center  of  gravity  and  elastic  center  of  the  TWS,  S$  =  tti(xg  —  %c)  and  Iq  denote 
its  static  and  polar  moments  of  inertia  about  the  elastic  axis,  and  Kh  and  Kq  are 
its  bending  and  torsional  stiffness  coefficients.  From  an  aerodynamic  viewpoint, 
a  TWS  has  the  shape  of  the  airfoil  of  the  wing. 

In  order  to  design  a  TWS  for  the  F-16  Block  40,  we  start  from  a  detailed 
FE  structural  model  of  a  clean  right  wing  of  this  fighter  with  a  missile  and 
its  launching  system  at  the  wing  tip  (see  Fig.  5).  We  determine  the  set  of 
parameters  Ptws  by  requiring  that  the  sought-after  TWS  be  equivalent  to  the 
wing  of  the  F-16  Block  40  in  the  following  sense 

1.  It  reproduces  the  first  bending  and  torsion  frequencies  of  the  dry  wing 
which  are  4.76  Hz  and  7.43  Hz,  respectively. 

2.  It  reproduces  the  same  bending  and  torsion  modal  masses. 

3.  It  reproduces  the  same  vertical  displacement  at  the  leading  edge  of  the 
cross  section  located  at  68%  of  the  distance  from  the  root  to  the  tip  of  the 
wing,  for  both  bending  and  torsion  mode  shapes,  when  these  are  normal¬ 
ized  to  a  unit  rotation  of  the  cross  section. 

For  criterion  3,  we  choose  the  wing  cross-section  located  at  68%  of  the  distance 
from  the  root  to  the  tip  (see  Fig.  5)  because  the  F-16  wing  is  strongly  tapered 
and  is  rather  soft  towards  its  tip.  The  criteria  stated  above  lead  to 


m 

=  2.05  103  Kg 

h 

=  2.53  103  Kg.m2 

xa 

=  1.126  m  (upstream  of  section  at  68%) 

(xG  -  xc) 

=  0.0642  m  (G  behind  C) 

Kh 

=  2.046  10®  N/m 

Ke 

=  5.468  10®  N.m 

We  choose  the  shape  of  the  typical  wing  section  as  that  of  the  airfoil  located 
at  45%  of  the  distance  between  the  root  and  tip  of  the  F-16  wing  (see  Fig.  5). 
This  choice  is  justified  by  two  factors:  (a)  the  chord  of  the  airfoil  of  the  TWS 
must  be  close  to  the  ratio  of  the  wetted  area  and  the  wing  span,  and  (b)  because 
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of  tapering,  most  of  the  lift  is  generated  by  the  section  of  the  wing  that  is  close 
to  the  root,  which  implies  that  the  aerodynamic  center  of  the  wing  is  within 
this  area.  We  discretize  the  two-dimensional  flow  domain  by  18,000  vertices, 
and  ensure  a  sufficient  resolution  for  shock  capturing  in  the  region  close  to  the 
thin  leading  edge  (see  Fig.  6). 

5.2  Detailed  FE  model 

Based  on  modeling  information  provided  by  Lockheed-Martin,  we  construct  a 
detailed  three-dimensional  FE  structural  dynamics  model  of  the  F-16  Block  40 
in  clean  wing  configuration  with  a  missile  and  launching  system  at  each  wing 
tip.  This  FE  model  features  bar,  beam,  solid,  plate,  shell,  metallic  as  well  as 
composite  elements,  and  a  total  of  168,799  dofs  (see  Fig.  7(a)).  Like  the  TWS, 
this  FE  model  reproduces  correctly  the  first  dry  bending  and  torsion  frequencies 
measured  to  be  4.76  Hz  and  7.43  Hz,  respectively.  Using  the  Gridgen™  software 
and  F-16  CAD  data  provided  by  the  Air  Force  Research  Laboratory  at  Wright 
Patterson,  we  ignore  the  wing  tip  missiles  and  generate  first  a  surface  grid  with 
63,044  grid  points  (see  Fig.  7(b)).  Then,  we  generate  a  fluid  volume  mesh  with 
403,919  vertices. 

5.3  Parameter  identification 

Let  fben  {ftor)  and  a&en  (ow)  denote  respectively  the  frequency  (in  Hz)  and 
damping  coefficient  of  the  first  aeroelastic  mode  of  the  F-16  configuration  consid¬ 
ered  herein  that  is  dominated  by  bending  (torsion).  To  determine  these  aeroe¬ 
lastic  parameters,  we  mimic  the  procedure  employed  in  flight  testing.  Whether 
using  the  TWS  or  detailed  FE  model,  we  excite  the  structure  in  an  appro¬ 
priate  manner  and  simulate  numerically  its  response  to  the  prescribed  initial 
disturbance.  In  the  case  of  the  TWS  model,  the  numerical  simulation  generates 
two  signals:  h(t)  and  8(t).  In  the  case  of  the  detailed  FE  model,  it  generates 
168,799  signals,  one  for  each  dof  of  this  model.  Hence,  in  the  latter  case,  we 
focus  only  on  the  vertical  displacement  dofs  at  one  tip  node  and  one  root  node 
of  each  wing.  This  corresponds  to  positioning  output  sensors  at  these  locations 
for  flight  testing.  Once  the  signals  are  generated  for  a  sufficiently  long  period 
of  time,  we  apply  the  Eigensystem  Realization  Algorithm  (ERA)  [36]  to  extract 
from  them  the  aeroelastic  parameters  ft,en ,  ftor ,  &ben,  and  a*or. 

The  ERA  is  a  real-time  parameter  identification  method  that  is  less  sensitive 
to  noise  than  the  classical  logarithmic  decay  method.  It  can  handle  multi-input 
multi-output  (multi-dof)  systems.  Most  importantly,  it  is  capable  of  identifying 
the  frequencies  and  damping  coefficients  of  the  two  lowest  dominating  modes  of 
a  structure  using  as  few  as  2  cycles  of  response  history,  as  long  as  the  sampling 
rate  is  on  the  order  of  500  to  1000  Hz  (A t  =  1  to  2  milliseconds)  [37].  Hence,  it 
is  particularly  attractive  to  time-domain  aeroelastic  applications  as  it  requires 
the  simulation  of  only  two  cycles  of  the  vibration  of  the  structure. 
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5.4  Validation 

We  have  chosen  the  F-16  configuration  described  in  this  paper  because  it  is 
the  only  configuration  for  which  we  were  able  to  obtain  modeling  information. 
However,  the  only  F-16  flight  test  data  we  were  successful  in  obtaining  is  for 
the  aeroelastic  torsional  mode  of  a  Block  40  configuration  flying  at  the  same 
altitude  of  3,000  m,  but  equipped  with  3  launchers,  2  pylons,  and  one  missile 
(close  to  the  tip)  under  each  wing.  Hence,  at  least  because  of  these  differences 
in  wing  configurations,  some  discrepancy  is  to  be  expected  between  the  results 
predated  by  our  simulations  for  0.8  <  Moo  <1-4  (for  example,  see  Fig.  8)  and 
the  available  flight  test  data. 

ihe  results  reported  in  Figs.  9-10  show  that 

•  The  aeroelastic  bending  frequencies  predicted  using  the  TWS  model  differ 
from  those  predicted  using  the  detailed  three-dimensioned  FE  model  by 
less  than  5%.  For  the  aeroelastic  torsional  mode,  the  maximum  relative 
difference  is  12.5%  and  occurs  around  M^  =  1. 

•  The  aeroelastic  torsional  frequencies  predicted  using  both  structural  mod¬ 
els  correlate  reasonably  well  with  flight  test  results.  Those  predicted  using 
the  detailed  three-dimensional  FE  model  are  however  more  accurate:  de¬ 
pending  on  the  Mach  number,  they  either  agree  very  well  with  the  aeroe¬ 
lastic  torsional  frequencies  measured  during  flight  test  or  differ  at  most  by 
7%. 

•  The  aeroelastic  damping  coefficients  of  the  torsional  mode  computed  us¬ 

ing  the  TWS  model  do  not  correlate  well  with  flight  test  data.  They 
erroneously  predict  flutter  for  0.9  <  <  1.4. 

•  The  aeroelastic  damping  coefficients  of  the  torsional  mode  obtained  using 
the  detailed  three-dimensional  FE  model  correlate  reasonably  well  with 
flight  test  data.  In  particular,  they  reproduce  the  same  trend  of  variation 
with  the  Mach  number. 

It  fellows  that  while  the  TWS  model  offers  obvious  practical  and  computa¬ 
tional  advantages  that  make  it  popular  —  at  least  in  the  research  community 
—  £;>r  flutter  analysis,  it  does  not  appear  to  be  reliable  for  fighter  applications. 
This  is  partly  because  the  TWS  model  is  realistic  only  for  fairly  homogeneous 
wings  with  high  aspect  ratios  and  small  angles  of  sweep.  However,  the  results 
reported  herein  also  show  that  when  applied  to  three-dimensional  structural 
and  fluid  models,  our  nonlinear  aeroelastic  simulation  capability  predicts  well 
the  aeroelastic  parameters  of  the  chosen  F-16  configuration,  including  in  the 
transonic  regime. 

5.5  Performance  results 

Finally,  we  briefly  discuss  the  speed  of  the  nonlinear  aeroelastic  simulation  tech¬ 
nology  described  in  this  paper,  particularly  in  the  context  of  this  second  quote 
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from  [6]:  “ Even  at  present,  existing  CFD  codes  should  be  able  to  obtain  five  flut¬ 
ter  solutions  in  one  year.  ”  For  this  purpose,  we  focus  on  the  more  interesting 
case  of  the  detailed  FE  model  of  the  F-16  fighter  with  168,799  dofs. 

First,  we  remind  the  reader  that  we  need  to  predict  only  the  first  two  cycles 
of  the  response  of  the  structure  to  an  initial  disturbance  (see  Section  5.3).  For 
the  F-16  configuration  considered  herein,  our  numerical  algorithms  sustain  a 
coupling  time-step  on  the  order  of  1  millisecond.  This  time-step  corresponds 
to  sampling  the  period  of  the  first  dry  torsional  mode  of  this  fighter  with  134 
points.  It  also  turns  out  that  this  time-step  is  such  that  271  time-steps  are 
needed  to  simulate  the  first  twc  cycles  of  the  structural  response.  Second,  we 
note  that  we  assume  a  linear  behavior  of  the  structure  (see  Eq.  (9)),  which  is 
justified  for  flutter  clearance  applications  in  level  flight. 

We  report  in  Table  1  the  performance  results  obtained  on  an  Origin  2000 
computer  equipped  with  R10000  195  MHz  chips,  as  a  function  of  the  number 
of  processors  Nproc  allocated  to  the  simulation.  These  results  correspond  to  a 
single  Mach  number  point,  and  two  cycles  of  the  response  of  the  structure.  In 
all  cases,  we  assign  a  single  processor  to  the  structure  solver  because  it  is  less 
computationally  intensive  than  the  flow  and  mesh  motion  solvers.  We  observe 
that  on  average,  60%  of  the  total  CPU  time  is  elapsed  in  the  flow  solver,  38% 
in  the  mesh  motion  solver,  and  only  2%  in  the  structure  solver.  The  small 
percentage  of  the  total  CPU  time  consumed  by  the  structure  solver  is  due  to 
the  assumed  linear  nature  of  the  structural  problem.  (We  note  however  that  for 
maneuvering  applications  where  nonlinear  geometric  effects  must  be  accounted 
for,  the  percentage  of  the  total  CPU  time  elapsed  in  structural  computations  is 
on  the  order  of  15%). 

The  parallel  speed-up  and  parallel  efficiency  results  reported  in  Table  1  high¬ 
light  the  good  parallel  scalability  of  our  three-field  based  nonlinear  aeroelastic 
simulation  technology.  One  can  reasonably  argue  that  today,  most  aerospace 
engineers  have  access  to  a  6-processor  computational  platform.  From  the  results 
reported  in  Table  1,  we  conclude  that  using  such  a  computing  system,  the  tran¬ 
sonic  aeroelastic  parameters  of  a  full  fighter  configuration  can  be  extracted  in 
12.8  hours.  For  a  given  Mach  number,  finding  the  flutter  speed  usually  requires 
a  bracketing  procedure.  It  is  our  experience  that  such  a  bracketing  procedure 
typically  incurs  4  to  5  simulation that  are  similar  to  the  one  discussed  herein. 
Hence,  for  a  given  Mach  number,  it  is  also  our  experience  that  extracting  the 
flutter  speed  requires  between  51  and  64  hours  CPU  on  a  6-processor  computa¬ 
tional  platform.  Therefore,  using  the  AERO-F,  AERO-S,  and  MATCHER  suite 
of  codes  on  a  6-processor  configuration,  we  can  obtain  five  flutter  point  solutions 
for  a  fighter  in  the  transonic  regime  in  less  than  2  weeks.  Using  24  processors 
for  this  purpose  reduces  the  total  simulation  time  to  less  than  4  days.  These 
estimates,  which  do  not  include  the  time  needed  for  building  the  FE  structural 
model  and  generating  the  fluid  grid,  are  also  confirmed  by  our  experience  with 
the  F-16  Block  40  aircraft. 
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Conclusions 


High  performance  military  aircraft  are  usually  flutter  critical  in  the  transonic 
regime  where  the  linear  flow  theory  fails  to  predict  correctly  the  unsteady  aero¬ 
dynamic  forces  acting  on  an  aircraft.  Consequently,  flutter  testing  of  scaled 
models  in  transonic  wind  tunnels  is  always  used  to  generate  corrections  to  flut¬ 
ter  speeds  predicted  by  linear  methods.  Because  the  design  of  a  wind  tunnel 
flutter  model  and  the  analysis  of  the  corresponding  data  require  over  a  year’s 
time,  it  has  been  suggested  that  computational  fluid  dynamics  (CFD)  based 
nonlinear  aeroelastic  simulations  could  be  used  as  a  replacement  for  wind  tun¬ 
nel  testing,  if  they  prove  to  be  practical  —  that  is,  fast  enough  —  and  reliable  [6]. 
We  believe  that  the  nonlinear  aeroelastic  simulation  methodology  presented  in 
this  paper,  as  well  as  other  similar  capabilities  developed  elsewhere,  axe  today 
sufficiently  mature  to  take  on  this  challenge.  At  the  present  time,  they  may 
not  be  sufficiently  fast  to  be  used  as  a  design  tool.  However,  on  a  6-processor 
computing  platform,  our  nonlinear  aeroelastic  simulation  technology  is  capable 
of  computing  five  flutter  point  solutions  for  a  fighter  in  the  transonic  regime 
in  less  than  2  weeks.  On  a  128-processor  system,  it  can  complete  the  same 
job  in  less  than  a  day.  In  general,  nonlinear  solution  methodologies  such  as 
the  one  overviewed  in  this  paper  do  not  target  flutter  problems  in  the  subsonic 
and  supersonic  regimes  where  simpler  linear  methods  appear  to  satisfy  industry. 
They  are  rather  meant  to  address  those  nonlinear  aeroelastic  problems  which 
simply  cannot  be  tackled  by  the  linear  theory  of  aeroelasticity.  These  include, 
among  others,  flutter  at  transonic  speeds,  buffeting,  and  prediction  of  transient 
loads  and  stresses  during  agressive  maneuvering.  The  nonlinear  aeroelastic  sim¬ 
ulation  capability  described  in  this  paper  has  tremendously  benefited  from  the 
technical  interaction  with  the  Flight  Test  Center  at  the  Edwards  Air  Force  Base, 
and  its  application  to  F-16  problems.  Similar  interactions  with  industry  should 
accelerate  its  transformation  into  a  “production”  tool. 
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Figure  2:  CSS:  the  conventional  serial  staggered  procedure 


Figure  3:  ISS:  the  improved  serial  staggered  procedure 
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Figure  4:  A  typical  wing  section:  elastic  center  (C),  center  of  gravity  (G), 
fuselage  (F),  angle  of  attack  (a0),  free-stream  velocity  (Voo)?  a  typical  point  in 
the  flow  domain  (N) 
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1:  cross-section  chosen  for  identifying  the  structural  properties 
2:  cross-section  chosen  for  identifying  the  aerodynamic  properties 


Figure  5:  FE  model  of  the  right  wing  of  an  F-16  Block  40 
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Figure  6:  Discretization  of  the  flow  domain  around  an  F-16  airfoil 
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(a)  Detailed  FE  structural  model 


(b)  Fluid  surface  grid 

Figure  7:  CSM  and  CFD  models  for  an  F-16  Block  40 
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Figure  8:  Mach  contours  and  streamlines  for  an  F-16  Block  40  at  Af< 


Mach  number 


Figure  9:  F-16  Block  40:  aeroelastic  frequencies  at  an  altitude  of  3,000  m 


(a)  Bending  Mode 


(b)  Torsion  Mode 


Figure  10:  F-16  Block  40:  aeroelastic  damping  coefficients  at  an  altitude  of 
3,000  m 
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Nproc 

CPU 

Total 

CPU 

Fluid 

CPU 

Mesh 

CPU 

Structure 

Parallel 

Speed-up 

Parallel 

Efficiency 

i 

69.8  hrs 

62.4  % 

37.4  % 

0.2% 

1.0 

100% 

3 

24.6  hrs 

64.6  % 

34.8  % 

0.6% 

2.8 

95  % 

6 

12.8  hrs 

63.3  % 

35.4  % 

1.3% 

5.4 

91  % 

12 

5.9  hrs 

57.1  % 

40.1  % 

2.8% 

11.9 

99  % 

24 

3.3  hrs 

52.7  % 

42.7  % 

4.6  % 

20.9 

87  % 

Table  Is  F-16  flutter  clearance:  Euler  flow  model  with  403,919  grid  points,  FE 
structural  model  with  168,799  dofs  —  Performance  results  on  an  Origin  2000 
for  a  single  Mach  number  point  and  two  cycles  of  the  response  of  the  structure 
(271  coupling  time-steps) 
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